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INTRODUCTION  -  subject.  Traumatic  brain  injury  (TBI)  is  a  major  contributor  to  morbidity  and 
mortality  experienced  by  those  soldiers  subjected  to  improvised  explosive  devices  (IED)  as  well  as  in  military 
and  civilian  high  speed  collisions.  Traumatic  cerebral  vasospasm  (TCV)  is  a  major  contributor  to  morbidity  and 
mortality  experienced  by  those  TBI  patients.  Aggressive  neurosurgical  treatment  motivated  by  early  diagnosis 
appears  to  improve  the  clinical  outcome  for  these  patients.  Early  transcranial  Doppler  (TCD)  measurements  of 
blood  flow  speed  within  major  cerebral  arteries  produces  the  initial  diagnosis,  hence  motivates  the  rapid 
treatment  of  TCV.  However,  the  skill  necessary  to  deploy  TCD  limits  its  availability  relative  to  its  need.  This 
gap  in  patient  care  reduces  the  quality  of  care  and  potentially,  therefore,  the  quality  of  life  of  injured  soldiers. 
This  gap  also  defines  a  critical  need  for  robust,  easy  to  use  TCD  system,  applicable  to  these  patients. 

INTRODUCTION  -  purpose.  Early  work  by  the  PI  on  automating  transcranial  Doppler  has  been 
translated  by  PhysioSonics,  Inc,  (PSI)  a  Seattle-based  company  that  he  co-founded  into  a  working  automated 
transcranial  Doppler  (aTCD)  device  called  ‘Presto’  that  they  have  submitted  to  the  FDA  via  a  51  OK  for 
approval.  Their  system  is  based  upon  a  novel  and  proprietary  ultrasound  platform  along  with  novel  and 
proprietary  ‘search  and  lock  on’  algorithms  for  detecting  and  staying  focused  on  the  brightest  spot  in  ‘Doppler’ 
space  -  that  is,  within  ultrasound-derived  maps  of  blood  flow  speed  captured  by  their  device.  This  is 
reasonable  because  a  view  of  each  of  the  major  cerebral  arteries  relevant  to  vasospasm,  listed  below,  shows 
that  each  contribute  the  largest  values  of  blood  flow,  hence  the  largest  signal  in  Doppler  space.  PhysioSonics 
has  optimized  it  to  measure  blood  flow  speed  within  the  middle  cerebral  artery  (MCA)  -  sufficient  for  some 
assays  of  vasospasm,  and  demonstrated  its  utility  on  uninjured  volunteers.  However,  a  critical  additional  assay 
of  vasospasm  -  the  Lindegaard  ratio  -  requires  measurements  of  blood  flow  speed  within  the  extra-cranial 
segment  of  the  internal  carotid  artery  as  well  as  from  the  MCA.  Moreover,  while  their  system  already  works  on 
humans,  PhysioSonics  has  no  plans  to  optimize  the  device  for  TCV,  or  to  extract  information  from  the  ICA. 

INTRODUCTION  -  scope  of  the  research.  To  met  the  goal  of  this  proposal  we  sought  to  test  the 
following  hypothesis:  PhysioSonics’  automatic  TCD  (aTCD)  system,  with  modifications  to  its  hardware  and/or 
software  as  necessary,  will  match  measurements  of  cerebral  blood  flow  in  the  middle  cerebral  artery  (MCA) 
and  in  the  internal  carotid  artery  (ICA)  made  by  standard  TCD  on  patients  who  suffer  from  cerebral  vasospasm. 

Specific  Aim  #  1 :  To  test  the  ability  of  the  aTCD  device  to  measure  blood  flow  velocity  in  the  MCA  and 
ICA  of  patients  with  known  vasospasm,  as  compared  to  standard  TCD  run  by  a  trained  neurosonographer. 
Proposed  study:  Use  the  device  to  measure  blood  flow  in  the  MCA  and  ICA  of  civilian  vasospasm  patients  - 
typically  those  with  sub-arachnoid  hemorrhage  -  within  Harborview  Medical  Center  (HMC),  the  only  Level  One 
trauma  center  in  Seattle.  Compare  those  measurements  with  those  determined  during  routine  clinical  care  by 
the  resident  neurosonographer  with  their  standard  clinical  TCD  unit.  Determine  the  level  of  fidelity  of  the 
aTCD-derived  information  on  blood  flow  with  that  given  by  the  gold  standard,  operator-dependent  TCD  device. 

Specific  Aim  #2:  To  modify  as  necessary  the  TCD  array  and  associated  holder  for  use  on  the  ICA, 
and,  the  algorithms  required  to  extract  the  highest  blood  flow  velocities  in  the  MCA  and  ICA  to  guarantee  that 
the  aTCD  device  can  detect  the  highest  levels  of  vasospasm.  Progosed  study  #1:  Modify  as  necessary  the 
brace  that  ordinarily  holds  the  transducer  to  the  temple  (to  facilitate  measurement  of  blood-flow  from  within  the 
MCA  so  that  it  will  automatically  extract  blood  flow  information  from  the  ICA.  Progosed  study  #2:  Modify  as 
necessary  the  software  within  the  aTCD  device  to  optimize  its  ability  to  calculate  the  value  of  fastest  blood  flow 
within  the  MCA  and  the  ICA  for  vasospasm  patients.  Progosed  study  #3:  Modify  as  necessary  the  TCD 
transducer  or  its  controlling  software  so  as  to  optimize  its  ability  to  locate  the  point  of  fastest  blood  flow  within 
the  MCA  and  the  ICA  for  vasospasm  patients 

Impact.  Successful  completion  of  the  proposed  work  will  yield  a  device  that  could  readily  become 
available  to  the  military  with  its  data  interpreted  on  site  or  sent  or  via  telemedicine  for  interpretation  by  distant 
neurosurgical  experts. 
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BODY 


Specific  Aim  #1 :  To  test  the  ability  of  the  aTCD  device  to  measure  blood  flow  velocity  in  the  middle 
cerebral  artery  (MCA)  and  internal  carotid  artery  (ICA)  of  patients  with  known  vasospasm,  as  compared  to 
standard  TCD  run  by  a  trained  neurosonographer.  (Q1-Q10  out  of  12  quarters) 

Task  1 :  Get  approval  for  the  human  subjects  study  (Q1  -2).  This  study  will  involve  up  to  1 28  patients 
with  vasospasm,  with  multiple  studies  per  patient,  over  the  three  years  of  the  study. 

(1.1)  Done  in  the  first  year. 

Task  2:  Take  possession  of  the  aTCD  system  (Q3). 

(2.1)  Done  in  the  first  year. 

Task  3:  Modify  the  software  associated  the  aTCD  system  so  it  will  report  and  store  the  spectrogram 
derived  from  each  of  the  MCA  and  ICA,  including  the  Lindegaard  ratio.  (Q2-5). 

(3. 1)  Done  in  the  first  year. 

Milestone  #1:  Achieve  ability  to  start  human  study.  Done. 

Task  4:  Modify  as  necessary  the  commercial  headband  that  holds  the  TCD  so  that  it  can  deploy  on  the 
upper  portion  of  the  neck  of  volunteers.  (Q1-2) 

(4. 1)  Done  in  the  first  year. 

Task  5:  Collect  then  archive  aTCD  and  standard  TCD  data  from  the  MCA  and  ICA  of  128  patients, 
along  with  appropriate  supporting  medical  data.  For  a  subset  of  these  patients  we  will  also  collect  standard 
TCD  data  twice  in  a  day  in  order  to  determine  the  length  of  time  we  will  need  to  eventually  monitor  patients  with 
our  aTCD  (Q3-Q10)  UW. 

(5.1)  We  started  to  recruit  patients  with  civilian  vasospasm  in  the  5/1 2th  quarter  of  our  study  (October 
2012).  Table  1  documents  the  number  of  patients  identified,  consented,  and  studied  during  year  two  of  this 
research  effort  (October  2012  -  September  15th  2013). 


Human  Data 

Identified 

Consented 

Studied 

#  of  data  sets 

Total 

1Oct2012- 

15Sept2013 

65 

22 

16 

25 

8th  of  12  quarters 
(1  June-1 5  Sept 
2013) 

33 

9 

11 

14 

Note  that  we  loose  two  thirds  of  our  identified  patients  to:  lack  of  bone  windows,  lack  of  available  next  of  kin, 
refusal  to  participate  (a  big  problem  because  these  patients  are  typically  awake  but  uncomfortable),  lack  of 
cranial  bone  and  patient  death.. 
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Task  6:  Compare  spectrogram  features  found  within  aTCD  and  standard  TCD  data.  (Q5-Q10)  UW 
(6. 1)  In  Figures  1-4  we  compare  representative  blood-flow  data  in  the  MCA  of  three  vasospasm 
patients.  (Note  that  we  get  pdfs  from  the  commercial  system  whose  resolution  for  this  report  fall  short  of 
adequate  representation  when  reproduced  here.  We’ll  improve  this  for  the  next  report.)  These  examples 
demonstrate  that  we  can  collect  TCD  data  from  the  MCA  of  patients,  with  a  mix  of  success  in  their  comparison 
with  the  gold  standard  data,  in  ways  that  have  or  will  produce  improvements  in  the  aTCD  system.  Our 
system’s  data  falls  short  in  two  key  areas  -  the  search  algorithm  and  the  spectral  envelope  calculation  -  as 
anticipated  during  the  construction  of  this  grant  application.  Indeed,  our  tasks  as  originally  designed  below 
seek  to  address  these  issues.  Also,  here  we  show  the  standard  TCD  system’s  reports  of  mean  flow  speeds, 
e.g.  107  cm/s  at  a  depth  in  the  left  MCA  of  50  mm,  and  pulsatility  index  (PI),  a  relative  measure  of  the  tendency 
of  the  peak  systolic  blood  flow  speed  to  rise  above  the  mean  blood  flow  speed.  The  aTCD  system  also  reports 
these  values,  but  in  a  different  window.  We  will  develop  quantitative  comparisons  of  these  quantities  as  well  as 
of  peak  systolic  blood  flow  speed.  Also,  our  ICA  measurements  with  the  aTCD  system  compare  favorably  with 
that  of  the  clinical  system.  We  have  focused  on  those  measurements  that  motivate  changes  to  the  aTCD 
system  here  for  this  report. 


Figure  1.  Data  from  the  aTCD  machine  (A)  and  a 
standard  TCD  machine  (B),  where  the  aTCD  data 
under-represents  blood  flow  speeds  in  the  MCA. 

Here  in  (A)  we  see  peak  systolic  velocities  identified 
by  the  commercial  spectral  envelope  varying  between 
1 00-1 50  cm/s.  However,  in  (B)  the  figures  show  that 
the  standard  TCD  machine,  with  user  intervention, 
measures  peak  systolic  velocities  as  a  function  of 
distance  from  the  transducer  [at  50,  56  and  42  mm 
distance]  with  values  closer  to  180  cm/s.  The  aTCD 
machine  provides  data  suggesting  this  patient  is  in 
mild  vasospasm  while  the  clinical  machine 
demonstrates  that  this  patient  is  in  severe 
vasospasm.  Therefore,  for  this  patient,  the  flow- 
envelope  calculator  fails  to  capture  peak  flow  values. 
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Figure  2.  Data  from  the  aTCD  machine  (A)  and  a 
standard  TCD  machine  (B),  where  the  aTCD  data 
well  represents  the  blood  flow  speeds  in  the  MCA 
as  defined  by  the  standard  TCD  machine.  Here  in 
(A)  we  see  peak  systolic  velocities  identified  by  the 
commercial  spectral  envelope  (in  fucia)  140-150  cm/s, 
nominally  close  to  the  underlying  spectrogram  (in 
green).  This  compares  favorably  to  measurements  of 
peak  systolic  blood  flow  speeds  (B)  by  the  standard 
TCD  machine  as  a  function  of  depth  [at  40,  44,  50  and 
54  mm  distance  from  the  transducer].  These  peak 
speeds  start  at  a  value  of  120  cm/s  at  40  mm  distance 
from  the  transducer,  increasing  to  close  to  150  cm/s 
by  54  mm  distance.  Both  the  aTCD  machine  and 
clinical  machine  demonstrates  that  this  patient  is  in 
moderate  to  severe  vasospasm.  Therefore,  for  this 
patient,  the  aTCD  system  performs  plausibly  well 
compared  to  the  clinical  system. 


B 


7 


L-MCA 


L-MCA 


4  2 

76.6 

Mean 

1.19 

pi 

•44  2 


98.2 


Mean 


1.32 


L-MCA 


mmmm^ 

pi 


L-MCA  60  -44-2 


L-MCA  64  -44  2 


L-MCA  60  *44  -  2 


Figure  3.  Data  from  the  aTCD  machine  (A)  and  a 
standard  TCD  machine  (B),  where  the  aTCD  data 
collects  representative  blood  flow  speeds  in  the 
MCA,  but  with  a  poor  spectral  envelope 
calculation,  and  at  too  shallow  a  depth,  as 
compared  to  the  data  collected  by  the  standard 
TCD  machine.  Here  in  (A)  we  see  peak  systolic 
velocities  identified  by  the  commercial  spectral 
envelope  (in  fucia)  110  cm/s,  woefully  low  compared 
to  the  150  cm/s  readily  identify  in  the  underlying 
spectrogram  (in  green).  These  ‘green’  values  of  blood 
flow  speed  compare  favorably  only  to  measurements 
of  peak  systolic  blood  flow  speeds  (B)  by  the  standard 
TCD  machine  at  its  shallower  depths  (between  44  and 
50  mm  distance  from  the  transducer  face).  The  aTCD 
search  algorithm  completely  misses  the  very  high  flow 
speeds  captured  by  the  standard  TCD  machine, 
namely,  at  or  over  240  cm/s  between  depths  of  50-64 
mm.  Therefore,  for  this  patient,  the  flow-envelope 
calculator  fails  to  capture  peak  flow  values  and  the 
device’s  search  algorithm  locks  onto  local  peak  flow 
velocities  that  are  too  shallow,  compared  to  the 
manually  steered  standard  TCD  system. 
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Task  7:  Write  up  results  for  publication  and  presentations.  (Q7-8;  Q10-12)  UW 

(7.1)  We  are  waiting  to  perform  this  task  until  we  implement  within  our  aTCD  system  optimal  flow- 
envelope  calculators  within  the  context  of  optimal  search  algorithms.  We  have,  however,  presented  these 
preliminary  results  to  our  neurosurgery  colleagues. 

Milestone  #2\  Develop  quantitative  comparisons  of  the  aTCD-derived  and  standard  TCD-derived  blood 
flow  data  from  vasospasm  patients,  thereby  highlighting  short  comings,  if  any,  of  the  aTCD  system  for  the  case 
of  vasospasm. 

We  have  collected  sufficient  number  of  data  sets  from  our  patients  to  identify  three  problems  with  our 
existing  aTCD  system:  [la]  the  flow-envelope  calculator  fails  to  capture  peak  flow  values  and  [1b]  the  device’s 
search  algorithm  locks  onto  local  peak  flow  velocities  at  too  shallow  a  depth  for  many  patients  with  vasospasm, 
all  relative  to  the  clinical  gold  standard.  In  addition,  thanks  to  the  formation  of  peripheral  edema  in  several  of 
these  patients,  we  can  say  that  [1c]  the  device’s  search  volume  is  too  shallow  for  an  appreciable  percentage  of 
patients  with  vasospasm.  Task  9  addresses  problem  #1  a  while  Task  10  addresses  problems  #1b,c. 

Milestone  #3:  get  feedback  on  the  quantitative  aspects  of  our  work  from  our  scientific  and  clinical  peers. 

Our  clinical  peers  are  quite  excited  by  the  prospect  of  successful  semi-automation  of  the  TCD  process 
but  recognize  that  we’ve  a  ways  to  go  before  we  succeed.  In  particular,  our  military  TCD  consultant,  Alex 
Razumovsky,  PhD,  has  verified  the  importance  of  this  approach  to  TCD  measurements  and  has  asked  us  to 
(a)  increase  the  depth  of  field  of  search  for  fastest  blood  flow,  to  (b)  improve  the  signal-to-noise  ratio  of  our 
system,  to  (c)  explore  use  of  an  intermediate  variable  -  a  time-varying  map  of  the  MCA  based  upon  creation  of 
an  ‘iso-surface’  of  Doppler-derived  blood-flow  speeds  -  a  ‘Doppler  map’,  discussed  below  -  to  map  the  blood¬ 
vessel  structure  sufficiently  to  visualize  the  MCA  in  vasospasm.  This  last  issue  is  particularly  important 
because  some  patients  after  blast  show  large  blood  flow  speeds  in  the  MCA  not  because  they  are  in 
vasospasm  but  because  their  hematocrit  is  sufficiently  low  that  the  resulting  reduced  effective  viscosity  of  the 
blood  allows  for  higher  than  normal  blood  flow  speeds. 

Specific  Aim  #2:  To  modify  as  necessary  the  TCD  holder  for  use  on  the  ICA,  and,  the  algorithms 
required  to  localize  as  well  as  extract  the  highest  blood  flow  velocities  in  the  MCA  and  ICA,  in  order  to 
guarantee  that  the  aTCD  device  can  detect  the  highest  levels  of  vasospasm.  (Q3-Q12  out  of  12  quarters) 

Task  8:  Modify  as  necessary  the  brace  that  ordinarily  holds  the  transducer  to  the  temple  (to  facilitate 
measurement  of  blood-flow  from  within  the  MCA)  so  that  it  will  also  automatically  extract  blood  flow  information 
from  the  ICA  of  patients.  (Q4-Q8)  PSI 

_ (8.  D  We  have  a  new  design  for  the  holder,  shown  below  in  Figure  4. _ 


knob  for  tightening 
headband. 


holder  of 

transducer  for  MCA 


Figure  4  shows  the  new  headset  design  that 
facilitates  simultaneous  collection  of  ipsilateral 
MCA  and  ICA  blood  flow  speeds.  Here  we  show 
a  headset  modified  from  commercial  TCD  headsets 
to  reduce  weight  and  reduce  the  complexity 
associated  with  setting  up  the  headset  while  also 
facilitating  measurements  of  ICA  blood  flow  speeds. 
The  slide  for  holding  the  MCA  transducer 
accommodates  not  only  the  standard  TCD 
placement  (anterior  part  of  temporal  window)  but 
also  the  minority  of  patients  whose  optimal  temporal 
window  lies  near,  above,  or  even  posterior  to  the 
ear.  The  articulated  arm  holds  the  ICA  transducer, 
with  a  mechanism  that  (a)  locally  articulates  as  well 
as  (b)  swivels  about  its  axis. 


Milestone  #4:  Deliver  device  that  holds  the  transducer  that  supports  the  measurement  of  blood  flow 
velocity  in  the  MCA  and  ipsilateral  ICA  for  its  deployment  in  Aim  #1.  Done.  We’ve  encountered  an  important 
practical  problem  for  this  study,  however,  namely  that  no  test  subjects  wish  to  use  the  headset,  given  the 
general  discomfort  of  our  patients,  all  of  whom  have  not  had  sufficient  subarachnoid  hemorrhage  to  require 
induction  of  a  therapeutic  coma.  The  headset  is  comfortable  enough,  as  tested  by  engineers.  The  vasospasm 
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patients  are  profoundly  uncomfortable,  however,  and  not  obliged  to  wear  the  headset  per  our  existing  IRB 
protocol.  We  will  talk  to  the  program  manager  about  adding  a  new  Task  to  this  study  where  we  will  refine  the 
designs  of  this  headset  using  feedback  from  healthy  volunteers. 

Task  9:  Modify  as  necessary  the  software  within  the  aTCD  device  to  optimize  its  ability  to  calculate  the 
value  of  fastest  blood  flow  (the  spectral  envelope)  within  the  MCA  and  the  ICA  for  vasospasm  patients.  (Q4- 
Q10)  UW  and  PSI 

First  recall  that  the  device  as  delivered  will  report  values  of  the  spectral  envelope.  Our  job  requires 
optimization  of  this  process  for  vasospasm  patients,  a  notoriously  difficult  cohort.  We  have  identified  two 
pathways  towards  development  of  optimized  means  of  calculating  the  value  of  fastest  blood  flow  in  the  MCA 
and  ICA,  referred  to  here  as  the  spectral  envelope:  a  statistical  analysis  and  a  mathematical  modeling  analysis. 
We  have  published  results  from  the  first  pathway  and  continue  to  refine  the  details  of  the  second  method. 

(9. 1)  Regarding  the  statistical  analysis  for  generation  of  plausible  spectrograms,  we  have  examined  a 
method  called  Gaussian  Mixture  Model  (GMM)  for  estimating  the  density  function  of  the  spectrogram  at  a  given 
time  slice.  Key  here,  we  are  trying  to  develop  an  independent  and  quantitative  metric  for  ‘maximum  blood  flow 
speed’  -  that  is,  for  the  spectral  envelope  used  by  the  clinicians  rather  than  ‘it  looks  good  by  eye’  which 
depends  heavily  on  the  choice  of  filter,  the  choices  of  color  for  representing  the  flow,  among  other  qualitative 
variables.  The  Journal  of  Ultrasound  in  Medicine  has  published  our  paper  on  this  subject  (Marzban  et  al,  2013; 
attached  at  the  end  of  this  report),  here  exercised  on  non-vasospasm  patients  for  which  there  exists  a  robust 
means  of  calculating  the  spectral  flow  envelope  using  a  commercial  algorithm  that  acts  as  our  gold  standard  for 
comparing  against  our  own  algorithm.  (Below  we  apply  our  algorithm  to  some  of  our  vasospasm  patients.) 

Here  we  summarize  our  results.  In  this  approach,  one  assumes  the  overall  distribution  of  blood  flow 
speeds  at  a  given  instant  in  time  consists  of  a  mixture  of  two  Gaussian  ‘density’  functions  where  density  refers 
to  the  blood-flow  speed  distribution  (Figure  5). 

Figure  5.  Top  panel.  Direct  measurements  of  the 
distribution  of  blood  flow  speed  within  the  middle 
cerebral  artery  (vertical  axis,  in  cm/s)  as  a  function  of 
time  (horizontal  axis),  here  for  four  beats  of  the  heart 
of  a  nortvasosgasm  patient  collected  during  another 
study.  The  various  lines  mark  the  instantaneous 
value  of  maximum  blood  flow  speed  throughout  the 
cardiac  cycle  generated  by  various  algorithms. 

Specifically,  this  spectrogram  of  a  non-vasospasm 
patient  (in  color),  shows  the  Spencer  envelope  (in 
black),  and  three  envelopes  according  to  the  GMM 
approach,  corresponding  to  the  90th,  95th,  and  99th 
percentiles  of  the  distribution  of  the  “signal"  defined  in 
the  bottom  panel.  In  that  bottom  panel  we  see  the 
spectrogram  (circles,  with  horizontal  axis  proportional 
to  flow  velocity  (the  vertical  axis  in  the  top  panel)  and 
the  vertical  axis  proportional  to  the  acoustic  power  in 
the  measurement  of  blood  flow)  at  one  time  slice  (at 
the  last  time  shown  in  the  top  panel),  with  specific 
definition  of  the  “background"  Gaussian  (red  curve), 
the  “signal"  Gaussian  (black  curve),  the  three 
percentiles  of  the  signal  (blue  vertical  lines),  and 
Spencer's  value  (red  vertical  line,  under  the  dark  blue 
line  associated  with  the  95th  percentile). _ 

The  specific  “mixing  ratio"  of  the  two  Gaussian  functions  and  the  individual  parameters  of  the 
Gaussians  are  estimated  from  data,  in  our  case  via  the  Expectation-Maximization  algorithm.  We  set  one  of  the 
two  Gaussian  components  to  interpolate  the  “background"  signal,  which  includes  the  noise  floor  of  the  signal 
as  well  as  occasional  transposition  of  retrograde  flow  to  very  high  (and  incorrectly  identified  as  such)  blood  flow 
speeds  as  in  Figure  5  above.  We  assign  the  second  Gaussian  component  with  the  actual  blood  flow  signal, 
from  its  low  values  near  the  blood-vessel  wall  to  its  maximum  speed  near  the  center  of  the  blood  vessel. 
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With  this  interpolation  of  the  distribution  of  blood-flow  speed  at  a  given  instant,  for  all  of  the  instants  of 
data  that  we  have  for  this  patient,  we  then  calculate  candidate  values  of  the  maximum  blood  flow  speed 
(shown  here  in  Figure  6). _ _ 


Figure  6.  Histograms  of  flow  velocity  (FV)  for  patient 
X  at  a  given  instant  in  time  (top),  and  the 
corresponding  cumulative  (ie,  integrated)  histogram 
(bottom).  Note  that  the  interpolative  Gaussian 
functions  are  not  shown.  When  the  cumulative 
frequency  reaches  a  pre-determined  value  (90,  95,  or 
99%  of  the  maximum  value),  we  define  the  flow 
velocity  there  as  that  associated  with  the  maximum 
flow  velocity.  In  the  bottom  figure,  we’ve 
demonstrated  this  process  for  the  95%  mark. 


_ fv _ | _ 

Based  upon  comparison  with  the  Spencer  flow-speed  envelope,  we  found  the  most  efficacious 
comparison  of  our  calculations  with  that  of  the  Spencer  algorithm  was  achieved  at  the  95th  percentile  of  the 
observed  blood  flow  speeds  by  our  GMM  approach  (Figure  7).  Specifically,  we  see  that  the  Spencer  envelope 
(that  is,  the  estimation  of  the  instantaneous  fastest  blood  flow  values  given  by  a  leading  commercial  algorithm) 
lies  nearest  to  the  95th  percentile  of  measured  blood  flow  speeds  calculated  by  our  algorithm. 
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As  reported  in  an  earlier  quarterly  report  we  have  applied  our  Double-Gaussian  spectral  flow  envelope 
calculator  to  a  subset  of  our  new  vasospasm  patients  (Figure  8)  with  success  thus  far,  even  when  we  introduce 
artificial  noise  into  the  data  as  a  means  of  testing  the  ability  of  our  blood  flow  spectral  envelope  calculator  to 
identify  the  fastest  blood  flow  under  conditions  consistent  with  severe  vasospasm  (low  signal-to-noise  in 

particular). _ 

Figure  8.  Here  we  demonstrate  the 
effect  of  artificially  added  noise  on  the 
ability  of  our  double-Gaussian  method 
(blue  lines)  to  calculate  the  flow 
envelope,  from  data  collected  from  a 
vasospasm  patient  using  our  Presto 
device.  Not  only  can  the  envelope 
eventually  fail  to  identify  the  point  of 
fastest  blood  flow  but  the  structure  of  its 
predictions  through  time  changes,  over¬ 
estimating  the  flow  at  a  greater  rate  at 
systole  than  throughout  the  rest  of  the 
time  series.  The  commercial  envelope 
(black)  consistently  under  predicts  the 
systolic  peak  in  blood  flow  speed,  a 
critical  failure  in  the  context  of  grading 
the  presence  and  extent  of  vasospasm. 

(9.2)  Regarding  the  mathematical  analysis  whose  foundation  is  direct  calculation  of  blood  flow  speeds 
in  compliant  blood  vessels  based  upon  extractable  information  from  the  measured  full  spectrogram,  that  is  the 
entirety  of  measurements  of  distributed  blood  flood  flow  in  the  vessel  through  time  from  a  given  patient.  The 
model  then  iterates  to  optimize  the  fit  between  the  mathematical  model  and  those  measurements.  Finally,  the 
given  the  mathematical  model  results  for  predicted  blood-flow  distribution  as  a  function  of  time  within  the  blood 
vessel,  a  noise-free  mathematically-based  estimation,  the  model  calculates  the  spectral  flow  envelope  for  that 
patient.  We  have  now  applied  it  to  some  of  the  human  data  we  collected  from  a  previous  study.  We  describe 
this  highly  detailed  mathematical  model  and  first  applications  to  human  data  in  a  draft  of  our  paper  on  the 
subject  (Morrison  et  al,  2013),  appended  to  this  report.  Here  we  summarize  its  most  salient  features  in  Figures 
9  and  10  below. 
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Figure  9.  Mathematical  model  for  predicting  blood-flow  parameters  from  noisy  measurements  of  the 
full  distribution  of  blood  flow  within  a  blood  vessel.  (A)  Summary  diagram  of  the  mathematical  model  and 
its  constituent  parts.  The  model  starts  with  the  observed  spectrogram,  that  is,  the  time  course  of  the  full 
distribution  of  blood  flow  within  the  blood  vessel  (from  zero  at  the  blood-vessel  walls  to  its  maximum  value 
within  the  center  of  the  blood  vessel).  Those  measurements  produce  first  plausible  values  of  initial  parameters 
that  drive  a  mathematical  model  for  blood-flow  distribution  (B)  in  compliant  blood  vessel  such  as  the  MCA  or 
ICA.  From  that  calculated  flow  the  algorithm  produces  first  the  spectrogram  associated  with  that  blood  flow 
(e.g.,  what  a  TCD  system  would  measure  from  this  predicted  blood  flow)  then  compares  that  result  with  the 
measured  spectrogram.  In  (C)  we  offer  three  different  intermediate  predictions  of  those  spectrograms  for  a 
given  patient,  each  based  upon  different  flow  parameters.  (Specifically,  in  C  we  show  the  observed 
spectrogram  in  the  top  row,  the  predicted  spectrogram  in  the  middle  row  and  the  associated  flow  field  in  the 
blood  vessel  in  the  bottom  row,  for  a  single  patient’s  data.)  After  iterating  in  this  fashion,  the  net  effect  is  a 
prediction  of  the  spectral  envelope  (among  other  medical  parameters  of  interest  such  as  pulsatility  index,  mean 
blood  flow  speed),  which  we  show  in  Figure  10,  below. _ 
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Figure  10.  (A)  Samples  of  summary  of  measurements  of  TCD-derived  blood-flow  distribution  (left  column), 
predictions  of  blood-flow  distribution,  and  predictions  of  TCD-derived  blood-flow  distribution  (right  column), 
from  5  different  patients.  We’ve  optimized  this  algorithm  to  match  to  the  best  of  its  abilities  the  blood-flow 
values  at  peak  values,  important  for  diagnostic  tests  in  the  case  of  vasospasm.  We  consider  these  results 
promising,  hence  anticipate  further  work  on  the  methodology;  nonetheless  this  algorithm  will  require  more  work 
before  its  deployment  in  the  clinic. 


Milestone  #5:  Deliver  software  that  robustly  measures  the  values  of  fastest  blood  flow  measured  within 
the  MCA  and  ICA  for  testing  and  validation  in  Aim  #1.  UW.  As  described  above  we  have  developed  two 
plausible  pathways  to  replace  existing  commercial  algorithms  for  calculating  the  fastest  blood  flow  with  new 
ones  optimized  for  vasospasm,  where  the  signal  to  noise  ratio  is  low  and  the  blood  flow  speeds  are  high,  each 
testing  the  limits  of  the  existing  system.  In  future  work  we  will  seek  additional  means  of  calculating  this 
quantity  as  well  as  test  these  new  methods  on  our  existing  data  sets. 
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Task  10:  Modify  as  necessary  the  TCD’s  controlling  software  so  as  to  optimize  its  ability  to  search  for 
and  locate  the  point  of  fastest  blood  flow  within  the  MCA  and  ICA  for  vasospasm  patients.  (Q4-Q1 0)  UW;  PSI 
(10.1)  As  reported  in  our  7th  quarter,  we  have  achieved  a  first  iteration  of  this  Task.  Our  studies  with 
patients  made  clear  the  need  to  have  the  ‘search  algorithm’  look  for  points  of  fastest  blood  flow  deeper  than  the 
standard  60  mm  within  each  of  the  middle  cerebral  artery  and  the  ICA.  We  have  made  that  change  (Figure  1 1 , 
reported  earlier)  and  reduced  the  noise  in  the  collected  data. _ 
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Figure  11.  This  figure  shows  the  modified  user 
interface  that  guides  the  user  to  setup  the  automatic 
TCD  to  locate  the  point  of  fastest  blood  flow  in  the 
middle  cerebral  artery.  Here  we  photo-document 
(with  a  smart-phone  camera)  that  the  device  can 
search  up  to  70  mm  deep  within  the  cranium.  The 
user  specifies  the  search  depth  (60  or  70  mm)  and 
then  lets  the  machine  collect  its  data. 


However,  as  reported  above,  despite  the  fact  that  our  algorithm  can  search  over  a  depth  range  up  to  70 
mm  within  the  cranium,  our  search  and  locate  algorithm  thus  far  tends  to  lock  on  to  plausible  blood-flow 
maxima  within  the  MCA  at  shallower  depths  (40-50  mm)  as  compared  to  the  depth  of  vasospasm  identified  by 
the  commercial  TCD  (50-65  mm),  which  facilitates  search  by  hand  for  the  global  maximum  in  blood  flow  speed 
along  the  MCA.  We  have  begun  to  work  with  PhysioSonics  to  figure  out  a  way  around  this  problem. 

(10.2)  We  started  to  determine  what  it  will  take  to  access  an  internal  variable  within  the  aTCD  system 
that  facilitates  identification  of  the  point  of  fastest  blood  flow  within  the  cranium  in  the  region  of  the  MCA,  by 
producing  a  time-varying  map  of  an  iso-surface  of  Doppler-derived  blood-flow  speed  (Figure  12).  The  potential 
benefits  the  projects  are  three-fold.  (1)  We  expect  to  refine  the  process  of  identification  of  fastest  blood  flow  by 
analyzing  how  this  system  collects  and  processes  this  internal  imaging  data.  (2)  We  expect  to  produce  the 
ability  to  visualize  vasospasm  by  directly  measuring  the  diameter  of  the  blood  vessels  through  analysis  of 
these  Doppler  maps.  (3)  Production  of  these  maps  may  allow  assessment  of  vasospasm  in  the  case  of  low 
hematocrit  (which  can  produce  large  values  of  blood  flow  speed  without  vasospasmO  by  superimposition  of  the 
identified  point  of  fastest  blood  flow  with  the  associated  vascular  anatomy. 

This  will  require  substantial  work  ‘under  the  hood’  to  make  these  data  readily  available  for  both  off-line 

and  real-time  analysis.  We  have  started  the  process  of  setting  these  internal  tasks. _ 

Figure  12.  Samples  of  instantaneous  three- 
dimensional  maps  of  blood  flow  in  the  cranium. 
Here  we  display  in  one  map  a  three-dimensional  view 
of  the  internal  carotid  artery  (in  red)  and  where  it 
meets  the  middle  cerebral  artery  (MCA,  in  blue). 

Each  color  displays  a  single,  absolute  value  of  blood 
flow  speed  towards  or  away  from  the  transducer, 
thereby  outlining  the  blood  vessel.  The  other  three 
images  offer  three  different  perspectives  on  the  MCA 
(in  red)  and  the  anterior  communicating  artery  (blue), 
as  well  as  locating  the  position  of  the  point  of  fastest 
blood  flow  towards  the  transducer  (green  dot). _ 


TCD-angio  identified! 

(may  be  able  to  skip  CT-  or  MRI-angio 


Milestone  #6:  Deliver  software  that  optimally  identifies  the  point  of  fastest  blood  flow  in  the  MCA  and 
ICA  for  vasospasm  patients  for  testing  and  validation  for  its  analysis  in  Aim  #1.  UW 

Task  1 1 :  Perform  pilot  study  involving  32  patients  that  tests  in  a  blinded  fashion  the  efficacy  of  the 
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aTCD  compared  with  a  standard  TCD  system.  (Q1 0-Q1 2).  UW 
(1)  NA 

Task  12:  Travel  to  a  military  medical  facility  under  the  sponsorship  of  Dr.  Armonda  to  determine  what 
changes,  if  any  need  to  be  made  to  the  FDA  approved  aTCD  system  for  optimal  deployment  to  military 
facilities.  Q10-12).  UW 
(1)  NA 

Task  13:  Write  up  results  in  a  final  report  for  the  sponsor,  for  publication,  and  within  a  follow-on 
proposal  that  will  field  trial  an  aTCD  device  at  a  military  hospital  involving  patients  with  blast-induced  TCV. 
(Q10-Q12).  UW 
(1)  NA 

Milestone  #7:  Report  results  to  the  sponsor  and  begin  the  transition  of  this  technology  to  the  military. 
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KEY  RESEARCH  ACCOMPLISHMENTS 


•  We  have  created  an  improved  headset  for  use  in  deploying  simultaneous  MCA  and  ICA  monitoring. 

•  We  have  verified  the  existence  of  three  technical  obstacles  between  the  current  system  and  one 
useful  for  vasospasm  detection  and  monitoring,  all  anticipated,  as  a  class,  in  writing  the  original 
proposal,  [la]  the  flow-envelope  calculator  can  fail  to  capture  peak  blood-flow  values;  [1b]  the 
device’s  search  algorithm  locks  onto  local  peak  flow  velocities  at  too  shallow  a  depth  for  many 
patients  with  vasospasm;  [1c]  the  device’s  search  volume  is  too  shallow  for  an  appreciable 
percentage  of  patients  with  vasospasm. 

•  To  address  problem  [1  a]  we  have  demonstrated  a  useful  means  of  calculating  spectral  envelopes 
based  upon  a  first-principles  approach  to  the  mathematical  modeling  of  blood  flow  in  the  blood 
vessel  of  interest. 

•  Also  to  address  problem  [1  a]  we  have  demonstrated  a  potentially  useful  means  of  calculating 
spectral  envelopes  based  upon  local  mathematical  information,  processed  in  a  statistical  way. 

•  To  address  [1c]  we  have,  earlier  this  year,  modified  the  internal  aTCD  software  in  order  to  increase 
the  depth  over  which  the  internal  search  algorithm  will  seek  to  identify  the  point  of  fastest  blood  flow 
in  the  MCA  and  ICA.  This  summer  we  verified  that  this  modification  works  while  in  use  on  patients. 

•  To  address  problem  [1b]  we  have  identified  an  approach  and  personnel  to  solve  this  problem. 

•  We  have  collected  25  data  sets  with  the  aTCD  system  from  1 6  patients,  far  short  of  our  yearly 
goals.  This  has  occurred  due  to  a  paucity  of  patients  to  consent  relative  to  historical  averages,  and 
because  only  one  third  of  the  available  patients  are  willing  to  participate  due  to  their  extreme 
discomfort.  To  address  this  problem  we  have  identified  a  parallel  site  within  the  Seattle  area 
(Swedish  Hospital  that  appears  willing  to  pursue  data  collection)  from  which  we  hope  to  collect  a 
comparable  number  of  data  sets  each  year. 

•  Also,  at  times  the  aTCD  system  requires  upgrades  that  take  it  out  of  commission.  To  continue  to 
collect  data  we  have  modified  our  IRB  (with  approval  now  from  the  local,  UW  committee)  to 
continue  to  collect  data  with  a  modified  standard  TCD  system,  which  will  allow  us  to  continue  to 
accrue  data  in  support  of  improving  the  spectral  envelope. 
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CONCLUSIONS 


•  We  can  collect  high  quality  TCD  data  using  the  aTCD  system  from  vasospasm  patients. 

•  In  doing  so  we  have,  however,  verified  the  existence  of  three  technical  obstacles  between  the 
current  system  and  one  useful  for  vasospasm  detection  and  monitoring,  all  anticipated,  as  a  class, 
in  writing  the  original  proposal,  [la]  The  flow-envelope  calculator  can  fail  to  capture  peak  blood-flow 
values;  [1b]  the  device’s  search  algorithm  locks  onto  local  peak  flow  velocities  at  too  shallow  a 
depth  for  many  patients  with  vasospasm;  [1c]  the  device’s  search  volume  is  too  shallow  for  an 
appreciable  percentage  of  patients  with  vasospasm. 

•  We  have  solutions  to  all  three  of  these  technical  problems,  in  various  states  of  implementation. 

•  We  have  also  identified  another  problem,  namely  a  reduced  number  of  patients  available  for  study, 
due  to  unusually  good  weather  in  Seattle  [civilian  subarachnoid  hemorrhage  often  occurs  in  Seattle 
during  times  of  bad  weather]  as  well  as  a  reluctance  of  these  patients  to  participate  in  the  exam 
because  they  are  awake  but  extremely  uncomfortable,  and  use  of  the  aTCD  (versus  a  standard 
TCD  system)  is  optional  from  a  clinical  perspective.  To  address  this  problem  we  have  identified  a 
parallel  site  within  the  Seattle  area  (Swedish  Hospital  that  appears  willing  to  pursue  data  collection) 
from  which  we  hope  to  collect  a  comparable  number  of  data  sets  each  year. 
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A  Double-Gaussian,  Percentile-Based 
Method  for  Estimating  Maximum 
Blood  Flow  Velocity 
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Objectives — Transcranial  Doppler  sonography  allows  for  the  estimation  of  blood  flow 
velocity,  whose  maximum  value,  especially  at  systole,  is  often  of  clinical  interest.  Given 
that  observed  values  of  flow  velocity  are  subject  to  noise,  a  useful  notion  of  “maximum” 
requires  a  criterion  for  separating  the  signal  from  the  noise.  All  commonly  used  criteria 
produce  a  point  estimate  (ie,  a  single  value)  of  maximum  flow  velocity  at  any  time  and 
therefore  convey  no  information  on  the  distribution  or  uncertainty  of  flow  velocity. 
This  limitation  has  clinical  consequences  especially  for  patients  in  vasospasm,  whose 
largest  flow  velocities  can  be  difficult  to  measure.  Therefore,  a  method  for  estimating 
flow  velocity  and  its  uncertainty  is  desirable. 

Methods — A  gaussian  mixture  model  is  used  to  separate  the  noise  from  the  signal  dis¬ 
tribution.  The  time  series  of  a  given  percentile  of  the  latter,  then,  provides  a  flow  veloc¬ 
ity  envelope.  This  means  of  estimating  the  flow  velocity  envelope  naturally  allows  for 
displaying  several  percentiles  (eg,  95th  and  99th),  thereby  conveying  uncertainty  in  the 
highest  flow  velocity. 

Results — Such  envelopes  were  computed  for  59  patients  and  were  shown  to  provide 
reasonable  and  useful  estimates  of  the  largest  flow  velocities  compared  to  a  standard 
algorithm.  Moreover,  we  found  that  the  commonly  used  envelope  was  generally  con¬ 
sistent  with  the  90th  percentile  of  the  signal  distribution  derived  via  the  gaussian  mix¬ 
ture  model. 
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Conclusions — Separating  the  observed  distribution  of  flow  velocity  into  a  noise  com¬ 
ponent  and  a  signal  component,  using  a  double-gaussian  mixture  model,  allows  for  the 
percentiles  of  the  latter  to  provide  meaningful  measures  of  the  largest  flow  velocities 
and  their  uncertainty. 

Key  Words — blood  flow;  brain;  head  injury;  noninvasive;  transcranial  Doppler  sonography 

Transcranially  derived  Doppler  measurements  of  blood  flow 
in  major  cerebral  arteries  have  found  many  clinical  applica¬ 
tions.1  In  addition  to  assaying  stroke  risk  due  to  sickle  cell 
disease,  dysfunction  of  cerebral  autoregulation,  and  a  patent  fora¬ 
men  ovale,  among  other  etiologies,  some  of  the  earliest  applica¬ 
tions  targeted  monitoring  for  vasospasm  after  subarachnoid 
hemorrhage.2-5  Cerebral  vasospasm,  the  transient  reduction  of 
the  diameter  of  1  or  more  major  cerebral  arteries,  can  lead  to 
reduced  blood  flow  into  the  brain  and  hence  cerebral  ischemia  and 
an  increasing  chance  for  neurologic  damage.1  Monitoring  for  the 
onset  of  vasospasm  remains  an  important  application  of  transcra¬ 
nial  Doppler  sonography,  with  more  30,000  patients  per  year  requir¬ 
ing  daily  monitoring  for  1  week  or  more.6  Adding  to  this  primarily 
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civilian  population  are  military  patients  with  head  injuries 
after  exposure  to  explosions,  typically  roadside  bombs, 
with  half  of  these  patients  having  cerebral  vasospasm.7 

Transcranial  Doppler  sonography  measures  the  dis¬ 
tribution  of  blood  flow  speeds  within  a  blood  vessel  toward 
or  away  from  the  transducer,  with  negligible  flow  speeds 
adjacent  to  the  blood  vessel  wall  and  maximum  flow 
speeds  near  the  center  of  the  vessel.8  Critical  for  the  assay 
of  cerebral  vasospasm,  among  other  uses  of  transcranial 
Doppler  sonography,  is  successful  capture  of  the  speed  of 
the  fastest  flowing  blood  within  the  major  cerebral  arter¬ 
ies,  since  this  value  acts  as  a  useful  proxy  for  blood  vessel 
narrowing.1  Capturing  this  speed  is  a  particularly  chal¬ 
lenging  problem,  since  vasospasm  reduces  the  volume  of 
blood  flow  while  accelerating  the  blood  flow  speed,9  hence 
reducing  the  target  for  transcranial  Doppler  measurements 
while  straining  against  the  upper  limits  of  ultrasound  data 
processing  due  to  aliasing.10 

The  time  series  of  a  flow  velocity  histogram  is  gener¬ 
ally  referred  to  as  a  spectrogram,  and  the  time  series  of  the 
maximum  flow  velocity  is  called  an  envelope.  Although  the 
spectrogram  conveys  a  great  deal  of  information,  the  enve¬ 
lope  is  often  the  only  quantity  a  clinician  examines.  This 
practice  is  reasonable  because  the  information  contained  in 
a  spectrogram  can  be  displayed  in  different  ways,  leading  to 
different  conclusions.  For  example,  Figure  1  shows  a  spec¬ 
trogram  with  2  different  color  schemes;  in  the  top  panel, 
the  spectrogram  is  linearly  related  to  the  color  scale, 
whereas  in  the  bottom  panel,  the  colors  are  proportional  to 
the  square  root  of  the  spectrogram. 

Figure  1.  Spectrograms  of  2  different  mappings  for  assigning  color  to 
the  flow  velocity  (FV)  values:  linear  (top)  and  sguare  root  (bottom). 


This  inherent  ambiguity  in  the  information  gleaned 
from  a  spectrogram  also  reflects  itself  in  the  corresponding 
envelope.  In  practice,  observed  flow  velocity  values  can 
range  from  1  to  300  cm/ s.  Therefore,  to  obtain  a  useful 
estimate  of  maximum  flow  velocity,  one  must  introduce 
some  criterion  that  defines  what  is  meant  by  maximum. 
Many  such  criteria  (standards  for  transcranial  Doppler 
measurements11-15)  are  based  on  the  cumulative  his¬ 
togram  of  flow  velocity.  Figure  2  shows  an  example  of  the 
relative  frequency  histogram  of  flow  velocity  (top)  and 
the  corresponding  cumulative  histogram  (bottom)  at  a 
specific  time  for  a  specific  patient  (hereafter,  patient  X); 
the  latter  is  obtained  by  integrating  the  former  from  left  to 
right.  One  may  define  the  maximum  flow  velocity  as  the 
value  at  which  the  cumulative  histogram  "levels  off,”  but 
there  exist  different  criteria  corresponding  to  different 
objective  measures  of  that  point.11-15  The  vertical  bar  in 
Figure  2  marks  the  flow  velocity  at  which  it  is  maximum 
according  to  the  modified  geometric  method.11'12  The  top 
panel  in  Figure  3  shows  the  spectrogram  and  the  modified 
geometric  method  envelope  (in  black)  for  patient  X. 

Given  that  the  proposed  algorithm  is  compared  with 
the  modified  geometric  method  algorithm,  a  brief  review 
of  the  latter  is  in  order.  The  modified  geometric  method 
algorithm  and  its  variants11'12  effectively  rotate  (clockwise 
and  about  the  origin)  the  cumulative  histogram  by  some 
amount.  The  effect  of  such  a  rotation  is  that  the  point  at 
which  the  cumulative  histogram  levels  off  translates  to  a 
point  at  which  the  rotated  cumulative  histogram  reaches 

Figure  2.  Histograms  of  flow  velocity  (FV)  for  patient  X  at  a  given  instant 
in  time  (top),  and  the  corresponding  cumulative  (ie,  integrated)  his¬ 
togram  (bottom). 
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the  maximum.  Given  certain  smoothness  constraints,  it  has 
been  shown  that  this  point  corresponds  to  the  maximum 
flow  velocity.12  Note  that  the  modified  geometric  method 
algorithm  generates  only  point  estimates  of  the  maximum 
flow  velocity  and  so  provides  no  natural  measure  of  uncer¬ 
tainty.  This  limitation  is  evident  in  the  cumulative  his¬ 
togram  in  Figure  2:  Although  the  modified  geometric 
method  provides  an  estimate  of  the  point  at  which  the 
cumulative  histogram  levels  off,  quantifying  the  corre¬ 
sponding  uncertainty  is  by  no  means  unambiguous.  Con¬ 
veying  uncertainty  is  important  because  different  levels  of 
uncertainty  call  for  different  clinical  actions. 

The  main  proposals  in  this  study  were  ( 1 )  to  develop 
a  new  means  of  separating  the  noise  and  the  signal  com¬ 
ponents  of  a  spectrogram  and  (2)  to  quantify  the  largest 
values  in  the  signal  using  the  concept  of  percentile.  The  nth 
percentile  of  a  quantity  x  is  defined  as  the  value  of  x  for 
which  n%  of  the  x  values  are  smaller.  For  example,  the  95  th 
percentile  offlow  velocity  is  the  flow  velocity  value  at  which 
95%  of  flow  velocity  values  are  smaller.  An  envelope  based 
on  the  95th  percentile  of  the  histogram,  then,  provides  a 
meaningful  estimate  of  the  “top  5%”  of  flow  velocities. 
Moreover,  2  (or  more)  envelopes  at  different  percentiles 
can  convey  some  notion  of  uncertainty.  In  this  study, 
envelopes  based  on  3  percentiles  were  considered:  90th, 
95th,  and  99th. 

Figure  3.  Top,  Spectrogram  for  patient  X  (color  background)  and  4  esti¬ 
mates  of  the  envelope:  modified  geometric  method  envelope  (black) 
and  90th,  95th,  and  99th  percentile-based  envelopes.  Bottom,  Modi¬ 
fied  geometric  method  envelope  (black),  95th  percentile-based  enve¬ 
lope  (red),  and  90th  and  99th  percentile-based  envelopes  (blue).  FV 
indicates  flow  velocity. 
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The  use  of  a  percentile  to  quantify  the  largest  flow 
velocity  presumes  that  the  histogram  is  the  correct  his¬ 
togram  offlow  velocity.  In  practice,  observed  spectrograms 
are  contaminated  by  various  types  of  noise.  One  of  these 
sources  of  noise  is  evident  in  the  small  “hump”  appearing 
on  the  right  side  of  the  (relative  frequency)  histogram  in 
Figure  2.  Flow  velocity  values  and  their  associated  per¬ 
centiles  are  useful  only  if  they  pertain  to  the  non-noise 
component  of  the  histogram,  hereafter  called  the  signal  dis¬ 
tribution.  In  fact,  the  90th  percentile  of  the  histogram 
shown  in  Figure  2  is  near  the  right-most  hump  and  there¬ 
fore  unrealistic.  Therefore,  to  compute  useful  percentile- 
based  envelopes,  one  must  first  disambiguate  the  signal 
and  noise  contributions  to  the  histogram.  To  that  end, 
here,  a  gaussian  mixture  model  with  2  components16,17  was 
used  to  represent  the  noise  and  signal  components,  respec¬ 
tively.  The  component  appearing  to  the  left  (closest  to  the 
origin)  was  defined  as  the  signal  distribution.  Armed  with 
the  signal  distribution,  the  3  percentiles  were  computed  at 
each  time,  and  their  time  series  was  computed  for  59 
patients.  The  3  percentile-based  envelopes  for  patient  X 
are  shown  in  Figure  3,  both  with  the  spectrogram  (top 
panel),  and  without  (bottom  panel). 

In  this  article,  the  details  of  this  percentile-based 
approach  for  estimating  an  envelope  are  presented.  It  was 
found  that  the  estimates  are  visually  consistent  with  the 
underlying  spectrograms,  and  the  modified  geometric 
method  envelope  is  approximately  consistent  with  the 
95th  percentile  envelope.  Also,  it  is  shown  that  percentile- 
based  envelopes  naturally  allow  for  displaying  envelope 
uncertainty. 

Materials  and  Methods 

Data 

Patient  data  for  this  preclinical  study  were  collected  from 
a  variety  of  hospitals  in  the  United  States,  following  a  study 
at  the  University  of  Washington  led  by  Dr  Mourad.  The 
data  included  arterial  blood  pressure  and  intracranial  pres¬ 
sure,  but  those  data  were  not  used  for  this  study;  only  sono- 
graphically  derived  transcranial  Doppler  spectra  were  used. 
Further  details  on  the  patient  data  maybe  found  in  an  arti¬ 
cle  by  Marzban  et  al.18  In  accordance  with  the  Institutional 
Review  Board  for  each  hospital,  informed  consent  was 
obtained  from  all  patients  or  their  families. 

Doppler  spectral  time  series  ofblood  flow  speed  in  the 
middle  cerebral  artery  (or  flow  velocity  -  flow  velocity) 
were  acquired  at  each  institution  via  clinically  approved 
transcranial  Doppler  units.  Data  acquisition  lengths  varied 
from  5  to  30  minutes.  All  retrospective  data  processing  and 
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analysis  were  conducted  at  the  Applied  Physics  Laboratory 
of  the  University  of  Washington.  Doppler  spectral  time 
series  data  collected  via  the  Applied  Physics  Laboratory's 
hospital  cohort  were  initially  digitized  at  125  Hz  in  time  and 
with  a  frequency  resolution  of  128  Hz  at  a  given  time. 
We  then  down-sampled  these  Doppler  time  series  to  a 
resolution  of  40  Hz  in  time,  sufficient  to  resolve  all  of  the 
following  features  in  all  of  our  patients'  data:  systolic  rise, 
diachrotic  notch,  and  diastolic  minimum.  We  also  selected 
a  fixed  duration  (118  time  steps,  or  ~3  minutes)  from  each 
of  the  59  patients  for  statistical  analysis.  This  duration  was 
sufficiently  long  to  include  several  cardiac  cycles  while  still 
allowing  details  of  the  envelopes  to  be  visually  evident. 


Data  Processing  and  Statistical  Analysis 
As  mentioned  in  the  introduction,  the  main  aims  of  the 
proposed  method  were  to  first  represent  the  instantaneous 
Doppler  frequency  distribution  of  the  flow  velocity  at  a 
given  moment  in  time  and  then  to  produce  spectral 
envelopes  of  those  frequency  distributions  through  time 
based  on  percentiles  of  that  distribution  in  the  Doppler  fre¬ 
quency.  Distinguishing  the  signal  from  the  noise  requires 
a  model.  In  this  article,  it  is  assumed  that  the  underlying 
distribution  offlow  velocity  consists  of  2  distributions,  cor¬ 
responding  to  the  signal  and  noise,  and  that  both  are  gauss- 
ian.  This  type  of  model  is  a  special  case  of  gaussian  mixture 
models,  wherein  a  distribution  is  assumed  to  be  a  linear 
combination  of  gaussian  distributions.16,17  The  weights 
in  the  linear  combination  (called  mixing  proportions)  and 
the  parameters  of  the  gaussian  distributions  are  then  esti¬ 
mated  from  data  via  some  optimization  procedure,  (in  this 
study,  the  expectation-maximization  algorithm  was  used, 
but  other  parameter  estimation  methods  are  equally  ade¬ 
quate.  The  expectation-maximization  algorithm  maxi¬ 
mizes  the  conditional  expected  log  likelihood.19)  Figure  4 
shows  the  signal  component  (black)  and  the  noise  com¬ 
ponent  (red)  for  the  example  shown  in  Figure  2.  Also 
shown  are  the  locations  of  the  maximum  flow  velocity 
according  to  modified  geometric  method  (black  vertical 
line)  and  the  90th,  95th,  and  99th  percentiles  of  the  signal 
distribution  (blue  vertical  lines).  All  statistical  analyses 
were  performed  in  R.20 

Results 


The  previous  section  illustrates  the  proposed  method  on  a 
given  patient  and  at  a  given  time.  The  top  panel  in  Figure 
3  shows  the  3  percentile-based  envelopes  (blue  curves)  for 
patient  X.  It  is  evident  that  the  percentile-based  envelopes 
are  consistent  with  the  underlying  spectrogram.  The  bot¬ 


tom  panel  in  Figure  3  shows  4  envelopes;  the  modified 
geometric  method  envelope  (black)  is  comparable  with  the 
95th  percentile  envelope  (red),  sandwiched  between 
the  90th  and  99th  percentile  envelopes  (blue).  Clearly, 
displaying  multiple  percentile-based  envelopes  conveys  a 
sense  of  the  uncertainty  in  the  maximum  flow  velocities. 

Although  technical,  it  is  also  worth  mentioning  2 
additional  steps  taken  to  generate  the  percentile-based 
envelopes  shown  in  Figure  3:  (l)  The  optimization  algo¬ 
rithm  used  for  the  gaussian  mixture  model  requires  initial 
estimates  of  the  parameters  of  interest.  At  time  0,  the  initial 
values  of  the  parameters  are  assigned  randomly,  and  the 
optimization  algorithm  produces  parameter  values  that 
characterize  the  2  best-fit  gaussian  components  of  the  flow 
velocity  histogram  at  that  time.  At  all  "future”  times,  the 
initial  values  of  the  parameters  are  set  to  the  final  values 
obtained  at  the  previous  time  step.  This  process  allows 
some  of  the  serial  (auto)  correlation  in  the  spectrogram  to 
be  reflected  in  the  envelopes.  If  the  values  of  the  gaussian 
mixture  model  parameters  are  initialized  randomly  for  every 
and  all  time  steps,  the  resulting  envelopes  are  somewhat 
less  smooth  than  those  shown  in  Figure  3.  (2)  A  running- 
median  filter  with  a  window  size  of  3  seconds  is  applied  to 
the  envelopes  to  smooth  them  even  further.  The  size  of  the 
window  (ie,  3)  is  obtained  by  trial  and  error  and  by  a  visual 
comparison  with  the  modified  geometric  method  enve¬ 
lope.  The  modified  geometric  method  envelope  too  is 
smoothed  by  a  variety  of  techniques,  eg,  an  averaging  (over 
multiple  cardiac  cycles)  filter  and  a  finite  impulse  response 


Figure  4.  Histogram  of  flow  velocity  (FV)  for  patient  X  (circles)  and  the 
estimated  densities  for  the  2  gaussian  components.  The  black  vertical 
bar  denotes  the  maximum  flow  velocity  according  to  the  modified  geo¬ 
metric  method,  and  the  blue  vertical  bars  correspond  to  the  90th,  95th, 
and  99th  percentiles  of  the  signal  distribution  on  the  left. 
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filter.  Although  the  smoothness  of  the  displayed  envelopes 
is  important  in  a  clinician  s  decision  making,  it  is  a  feature 
that  is  easily  controlled  (eg,  by  a  single  parameter,  such  as 
the  window  size  of  the  running-median  filter)  and  so  is  not 
of  serious  concern. 

The  procedure  was  applied  to  all  59  patients  in  the 
data  set.  It  is  impractical  to  show  all  59  figures  analogous  to 
Figure  3.  Instead,  the  3  percentile-based  envelopes  were 
compared  with  the  modified  geometric  method  envelope 
using  scatterplots,  which  were  then  further  summarized 
into  scalar  measures. 

Each  panel  in  Figure  5  shows  the  scatterplot  of  one  of 
the  percentile-based  envelopes  versus  the  modified  geo¬ 
metric  method  envelope  for  patient  X.  There  are  1 18  points 
on  each  scatterplot,  corresponding  to  the  118  time  steps 
displayed  in  Figure  3.  The  diagonal  line  is  simply  a  line  of 
slope  1,  y-intercept  0.  It  is  clear  from  the  middle  panel  that 
the  95th  percentile  envelope  generally  agrees  with  the  mod¬ 
ified  geometric  method  envelope.  As  expected,  the  90th  and 
99th  percentile  envelopes  are  generally  lower  and  higher, 
respectively,  than  the  modified  geometric  method  enve¬ 
lope.  The  disagreements  between  the  percentile-based 
envelopes  and  the  modified  geometric  method  envelope 
are  largest  for  higher  flow  velocity  values.  This  result  is  espe¬ 
cially  evident  in  the  99th  percentile  envelope  (bottom 
panel)  in  the  manner  in  which  the  deviation  of  the  scatter¬ 
plot  from  the  diagonal  line  is  largest  for  higher  flow  veloc¬ 
ity  values.  In  other  words,  the  systolic  peaks  of  the  99th 
percentile  envelope  are  in  fact  higher  than  might  be 
expected  from  the  modified  geometric  method  envelope. 

One  may  further  summarize  each  of  these  panels  by 
computing  the  root  mean  squared  error  between  the  mod¬ 
ified  geometric  method  envelope  and  each  percentile- 
based  envelope.  For  patient  X,  the  root  mean  squared  error 
values  corresponding  to  the  3  panels  in  Figure  5  are  17.9, 
16.9,  and  31.2  cm/s.  In  addition  to  confirming  that  the 
modified  geometric  method  envelope  is  closest  to  the  95th 
percentile  envelope  (for  patient  X),  these  numbers  have 
useful  interpretations  as  well;  eg,  it  can  be  said  that  the 
modified  geometric  method  envelope  and  the  95th  per¬ 
centile  envelope  have  a  typical  deviation  of  about  16.9 
cm/ s  across  the  118  time  steps  examined. 

Although  useful  for  comparing  different  envelopes, 
root  mean  squared  error  values  do  not  assess  whether  the 
difference  in  the  envelopes  is  due  to  the  amount  of  scatter 
in  the  scatterplot  or  to  an  overall  shift.  There  exists  a 
decomposition  of  root  mean  squared  error  that  allows 
one  to  examine  both  contributions.  The  details  of  the 
decomposition  are  presented  in  “Appendix:.”  Here,  suf¬ 
fice  it  to  say  that  the  correlation  coefficient  (denoted 


Corr)  and  bias  (ie,  mean  [modified  geometric  method]  - 
mean  [percentile-based  method] )  gauge  the  amount  of 
scatter  and  the  overall  shift,  respectively.  The  Corr  values 
for  the  3  panels  in  Figure  5  are  0.80,  0.81,  and  0.81  (from 
top  to  bottom),  suggesting  that  the  amount  of  scatter  in 
the  3  panels  is  comparable.  Said  differently,  any  of  the  3 
percentile-based  envelopes  provides  an  adequate  approx¬ 
imation  to  the  modified  geometric  method  envelope, 
if/when  correlation  is  the  measure  of  agreement.  The  main 
difference  between  the  3  panels  is  in  the  Bias  values:  8.5,  - 
2.6,  and  -23.5  (from  top  to  bottom).  These  numbers  have 
useful  interpretations,  too;  eg,  one  can  say  that  the  95th  per¬ 
centile  envelope  is  generally  shifted  above  the  modified  geo¬ 
metric  method  envelope  by  about  2.6  cm/ s  for  patient  X. 

To  compare  the  various  envelopes  across  all  of  the 
patients,  root  mean  squared  error,  Corr,  and  Bias  were 
computed  for  all  59  patients.  Their  distributions  are  shown 
in  Figure  6.  The  left-most  box  plot  in  the  top  panel  sum¬ 
marizes  the  histogram  (across  59  patients)  of  the  root 
mean  squared  error  computed  for  the  90th  percentile 
envelope  and  the  modified  geometric  method  envelope 
(across  118  time  steps). 

The  middle  panel  in  Figure  6  is  analogous  to  the  top 
panel  except  that  the  measure  of  agreement  between  the 
modified  geometric  method  envelope  and  the  percentile- 
based  envelope  is  Corr.  The  3  box  plots  are  comparable, 
implying  that  all  3  percentile-based  envelopes  are  compa¬ 
rably  correlated  with  the  modified  geometric  method 
envelope  across  all  59  patients.  The  bottom  panel  shows 
that,  on  average  (across  all  59  patients),  the  90th  and  99th 


Figure  5.  Scatterplots  of  the  3  percentile-based  envelopes  and  the 
modified  geometric  method  envelope  for  patient  X.  FV  indicates  flow 
velocity. 
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percentile  envelopes  are  below  and  above  the  modified 
geometric  method  envelope,  respectively.  By  comparing 
the  3  box  plots  in  the  top  or  bottom  panel,  it  can  be  seen 
that  the  agreement  between  the  modified  geometric 
method  envelope  and  the  95th  percentile  envelope,  noted 
for  patient  X,  extends  to  all  59  patients. 

It  is  worth  displaying  the  envelope  for  a  patient  for 
whom  such  scalar  measures  are  generally  poor.  Figure  7 
shows  the  spectrogram  and  envelopes  for  the  patient  with 
the  lowest  Corr  value.  It  is  clear  that  the  modified  geo¬ 
metric  method  envelope  (black)  is  in  fact  incorrect  and 
that  the  percentile-based  envelopes  (blue)  are  more  con¬ 
sistent  with  the  underlying  spectrogram.  In  other  words, 
the  low  Corr  value  for  this  patient  does  not  imply  a  defect 
on  the  part  of  the  percentile-based  envelope  but  rather  a 
defect  in  the  modified  geometric  method  envelope. 

Discussion 

It  has  been  shown  that  a  double-gaussian  mixture  model 
can  be  used  to  identify  a  component  of  the  histogram 
of  flow  velocity  corresponding  to  a  "signal.”  In  turn, 
percentiles  of  this  signal  distribution  provide  envelopes, 
which  meaningfully  quantify  the  largest  flow  velocity. 
Not  only  are  these  percentile-based  envelopes  visually 
consistent  with  their  underlying  spectrograms,  they 
are  also  objectively  consistent  with  a  commonly  used 
envelope  based  on  the  modified  geometric  method.  The 
percentile-based  envelopes  not  only  objectively  quantify 

Figure  6.  Distribution  of  the  root  mean  squared  error  (RMSE),  correla¬ 
tion  coefficient  (Corr),  and  Bias  (defined  in  “Appendix”)  for  all  59  patients 
examined  in  this  work. 


what  is  meant  by  "largest  flow  velocities”  but  also  have  the 
added  advantage  (over  existing  envelopes)  of  allowing 
for  a  natural  display  of  uncertainty  in  the  envelopes. 

Given  that  the  analysis  is  performed  on  real  (not  sim¬ 
ulated)  data,  the  true  envelope  is  not  known.  As  such,  the 
quality  of  each  envelope  cannot  be  assessed  objectively  on 
its  own  merit.  For  that  reason,  the  analysis  here  has  been 
based  on  visually  comparing  an  envelope  with  the  under¬ 
lying  spectrogram  (eg,  as  in  Figure  3),  or  objectively  com¬ 
paring  each  of  the  percentile-based  envelopes  with  the 
modified  geometric  method  envelope  (eg,  Figures  5  and 
6).  In  spite  of  the  central  role  played  by  the  modified 
geometric  method  envelope  in  comparing  envelopes,  it  is 
important  to  recall  that  the  modified  geometric  method 
envelope  is  not  the  true  envelope  (eg,  see  Figure  7).  Here, 
it  is  used  as  a  standard  only  because  of  its  common  usage. 

Of  the  3  percentile-based  envelopes,  the  95th  per¬ 
centile  envelope  is  closest  to  the  modified  geometric 
method  envelope,  but  the  agreement  is  not  perfect,  as  seen 
by  the  box  plots  in  Figure  6.  For  example,  according  to  the 
middle  box  plot  in  the  bottom  panel,  on  average  (across 
patients),  the  95th  percentile  envelope  is  slightly  above  the 
modified  geometric  method  envelope.  One  may  ask  what 
is  the  exact  percentile  corresponding  to  the  modified  geo¬ 
metric  method  envelope,  but  the  question  itself  is  inap¬ 
propriate  because  it  elevates  the  status  of  the  modified 
geometric  method  envelope  to  a  "reference  standard,” 
when  in  fact  it  is  not.  What  is  more  important  is  that  the 
proposed  approach  produces  a  useful  spectrogram  as  well 

Figure  7.  Same  as  Figure  3,  but  for  a  patient  for  whom  the  modified  geo¬ 
metric  method  envelope  (black)  is  clearly  incorrect.  FV  indicates  flow 
velocity. 
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as  an  interpretable  (in  terms  of  percentiles)  generalization 
of  approaches  that  naturally  produce  only  a  single  enve¬ 
lope  (eg;  modified  geometric  method). 

A  percentile-based  envelope  is  based  on  a  percentile  of 
the  signal  distribution;  in  which  the  signal  (and  noise)  dis¬ 
tribution  is  defined  from  fitting  a  gaussian  mixture  model 
with  2  components  to  the  whole  histogram;  at  a  given  time. 
There  is  some  justification  for  examining  the  3  compo¬ 
nents.  For  example;  in  Figure  4,  one  may  argue  that  there 
are  possibly  3  overlaying  histograms:  1  corresponding  to 
the  large  hump  on  the  left  (ie;  the  signal);  1  associated  with 
the  small  hump  on  the  right  (ie;  associated  with  an  alias¬ 
ing  reflection  on  Doppler  imaging;  and  another  corre¬ 
sponding  to  a  uniform  "background”  spanning  the  full 
range  of  flow  velocity  values.  We  repeated  the  entire  analy¬ 
sis  but  with  3  gaussian  components.  The  results  are  incon¬ 
clusive  and  require  further  research.  For  some  patients;  the 
results  do  not  change  substantially;  but  for  others;  they  do 
in  ways  that  can  be  considered  "better”  or  "worse”  depend¬ 
ing  on  the  measure  of  quality. 

The  percentiles  of  the  signal  distribution  allow  for 
meaningful  envelopes;  which  can  objectively  assess  the 
largest  flow  velocity.  In  addition;  displaying  multiple 
envelopes  can  convey  information  on  the  uncertainty  in  the 
observed  flow  velocity.  One  can  envisage  an  alternative 
measure  of  uncertainty.  For  example;  one  may  consider  the 
envelope  corresponding  to  a  fixed  percentile;  say  95th;  and 
then  obtain  the  distribution  of  that  quantity  via  some 
resampling  method.21  In  turn;  that  (sampling)  distribution 
can  be  used  to  compute  confidence  intervals  for  the  true 
95th  percentile  envelope.  Such  interval  estimates  of  the 
envelope  can  be  useful  in  conveying  uncertainty  if/ when  a 
specific  percentile  is  of  interest.  Otherwise;  displaying  mul¬ 
tiple  envelopes  corresponding  to  different  percentiles;  as 
done  here;  is  sufficient  for  conveying  uncertainty. 

Appendix 

Let  x.  and  y.  denote  the  values  of  the  modified  geometric 
method  envelope  and  a  percentile-based  envelope;  respec¬ 
tively;  at  time  i.  The  root  mean  square  error  between  2 
envelopes  is  defined  as  the  square  root  of  the  mean  square 
error  ( MSE ) : 

(1)  mse  =  T)  (  )2 

n  — 1  v  1 

where  n  is  the  length  of  the  time  steps  (here;  1 18).  It  can  be 
shown  that 

(2)  MSE  =  s2  +  s2  -  2(Corr \sy  +  /L  (Bias)2, 


where  s  and  s  are  the  sample  standard  deviations  of  x  and 

x  y  -t 

y,  respectively;  Corr  denotes  the  Pearson  correlation  coef¬ 
ficient;  and 

(3)  Bias  =  mean  x  -  meany. 

According  to  Equation  a  comparison  of  2  envelopes  in 
terms  of  mean  square  error  (or  root  mean  square  error)  is 
equivalent  to  a  comparison  in  terms  of  Corr  and  in  terms 
of  Bias;  the  former  compares  the  correlation  between  the 
envelopes;  and  the  latter  measures  the  overall  shift  between 
them. 
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Abstract 

A  physiological  model  of  blood  flow  in  the  middle  cerebral  artery 
is  put  forth.  The  model  consists  of  a  simple  circulatory  system,  pul¬ 
satile  flow  in  an  elastic  tube,  and  the  conversion  of  flow  velocity  to 
spectrograms  which  can  in  turn  be  compared  to  data.  Many  of  these 
component  models  have  parameters  which  are  estimated  by  measure¬ 
ments  of  blood  flow  with  transcranial  Doppler  ultrasound.  In  addition 
to  producing  spectrograms  which  closely  resemble  observed  spectro¬ 
grams,  the  model  provides  a  physical  basis  for  the  interpretation  of 
spectrograms. 


1  Introduction 

Modeling  blood  flow  has  been  a  topic  of  extensive  research  (refs??).  It  is 
important  not  only  for  the  purpose  of  understanding  the  underlying  dynamics 
(refs??),  but  also  for  its  diagnostic  value  (refs??).  Pierre?? 

Transcranial  Doppler  Ultrasound  (TDU)  provides  one  method  for  mea¬ 
suring  blood  flow  (refs??).  The  measurements  are  generally  displayed  in 
the  form  of  a  spectrogram.  A  spectrogram  is  a  measure  of  the  intensity  of 
acoustic  return  to  TDU  as  a  function  of  time  and  velocity,  where  velocity  is 
derived  from  signal  processing  to  extract  the  return  due  to  a  velocity  toward 
the  transducer.  A  spectrogram  is  a  function  of  time  in  the  sense  that  it  is 
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recorded  over  time.  Such  spectrograms,  however,  are  generally  corrupted 
by  a  number  of  errors,  including  measurement  errors,  return  from  unwanted 
sources,  etc.  Pierre??  A  basic  assumption  in  examining  Flow  Velocity  (FV) 
data  is  that  there  exists  an  underlying  spectrogram  which  adheres  to  certain 
physical  constraints.  For  example,  at  a  given  time,  the  brightness/amplitude 
of  a  spectrogram  is  expected  to  drop  to  zero  at  a  well-defined  (maximum) 
FV  value;  but  observed  spectrograms,  usually  do  not  reflect  that  physical 
requirement.  Also,  observed  spectrograms  are  often  contaminated  by  ran¬ 
dom  noise,  i.e.,  random  errors  which  cannot  be  attributed  to  any  physical 
constraints.  A  physics-  (or  physological)-based  model  is  required  to  allow 
one  to  attribute  physical  meaning  to  an  observed  spectrogram. 

The  approach  proposed  in  this  paper  involves  both  statistical  and  phys¬ 
ical  models.  Statistical  models  usually  involve  parametric,  mathematical 
expressions  whose  structure  is  dictated  by  the  convenience  of  estimating 
their  parameters.  Regression  models  are  canonical  examples  of  such  mod¬ 
els  (refs  ??;  Marzban  et  al.  2013).  By  contrast,  physical  models  begin  with 
physiologically-relevant  components  whose  mathematical  properties  are  well- 
understood,  and  attempt  to  describe  the  overall  phenomenon  (e.g.,  blood 
flow)  by  an  assembly  of  such  components. 

Physiologically  based  blood  flow  models  of  the  type  examined  here  have 
been  considered  in  the  past,  but  the  model  parameters  have  not  been  esti¬ 
mated  from  an  entire  spectrogram  (refs??).  For  example,  [15]  et  al.  model 
blood  flow  through  a  constriction,  but  at  no  stage  in  that  work  a  spectrogram 
is  used  to  estimate  the  model  parameters.  By  contrast,  when  a  spectrogram  is 
fully  employed,  the  focus  of  the  work  is  an  envelope,  and  not  a  2d  flow  veloc¬ 
ity  field  (refs??  geo  method,  double-gaussian  ...??).  One  novel  feature  of  the 
proposed  method  is  the  manner  in  which  a  full  physiologically-based  blood 
flow  model  is  utilized,  and  its  parameters  are  estimated  from  an  entire  spec¬ 
trogram.  The  estimation  method,  a  simple  instance  of  Recursive  Bayesian 
Estimation  (RBE),  (refs??)  involves  iteratively  making  observations,  and  re¬ 
vising  estimates  about  the  system  from  which  those  observations  originate. 
RBE  has  been  proven  useful  in  a  wide  range  of  applications  (refs??,  include 
subs),  including  the  modeling  of  blood  flow  [4]  (more  refs??). 

Note  to  self:  Cite  Windkessel  model  as  in  [2].?? 
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The  flow  chart  in  Figure  1  highlights  the  major  components  of  the  method¬ 
ology. 

First,  an  observed  spectrogram  is  consulted  in  order  to  choose  initial  values 
for  some  parameters.  These  values  are  then  used  in  a  model  of  the  circulatory 
system,  which  in  turn  produces  a  one  dimensional  time  series  of  maximum 
FV.  Solutions  to  approximations  of  Navier  Stokes  equations  in  an  elastic  tube 
are  used  to  convert  the  1-D  time  series  into  a  2-D  flow  field.  Although  the 
modeled  axial  and  radial  flow  velocity  and  vessel  wall  motion  can  be  used  for 
multiple  purposes,  the  main  attribute  of  interest  is  the  (ideal)  spectrogram 
computed  from  it.  All  of  the  component  models  which  lead  to  this  ideal 
spectrogram  are  referred  to  as  the  forward  model.  The  ideal  spectrogram  is 
then  compared  to  the  observed  spectrogram,  the  result  of  which  is  a  posterior 
Probability  Density  Function  (PDF)  of  some  of  the  parameters,  given  an 
observed  spectrogram.  From  this  PDF,  optimum  values  for  some  of  the 
parameters  are  obtained.  From  the  entire  model  a  number  of  clinically  useful 
quantities  can  be  generated,  such  as  envelope,  mean  flow,  or  pulsatility  index. 
The  complete  process  of  estimating  the  dynamic  and  fixed  parameters  is  often 
referred  to  as  inverting  the  forward  model. 


2  Data 


Ninety  six  recorded  spectrograms  are  examined,  73  of  which  are  collected 
with  the  Spencer  Technologies  TDU  device  in  a  hospital  setting.  The  major¬ 
ity  of  these  patients  have  experienced  traumatic  brain  injury,  and  five  have 
exhibited  vasospasm.  The  spectrograms  for  the  remaining  23  patients  are 
from  a  Presto  TDU  device  currently  under  development.  Pierre?? 

The  spectrograms  for  the  96  patients  are  complex  and  varied  across  pa¬ 
tients.  Although  the  spectrogram  for  each  patient  appears  to  be  consistent 
with  what  one  might  expect  from  an  ideal  spectrogram,  no  simple  ideal  spec¬ 
trogram  emerges  as  a  viable  candidate  when  multiple  patients  are  examined. 
In  fact,  the  spectrograms  for  9  of  the  73  patients  are  sufficiently  complex  to 
have  failed  a  commonly-used  envelope-detection  algorithm,  referred  to  as  the 
modified  geomteric  algorithm  [12]. 

Figure  2  shows  the  spectrogram  as  measured  by  a  TDU  device  (spencer  or 
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Figure  1:  In  the  proposed  methodology  the  forward  model  consists  of  several 
components,  each  with  parameters  which  are  to  be  estimated  from  observed 
spectrograms.  The  parameters  are  randomly  initialized,  and  the  resulting 
output  of  the  forward  model  is  compared  with  an  observed  spectrogram.  An 
RBE  approach  is  then  employed  to  update  the  parameters,  with  the  goal  of 
improving  the  agreement  between  observed  and  predicted  spectrograms. 


presto??).  Note  the  varying  shading  in  the  “body”  (i.e.,  the  lower  portion)  of 
the  spectrogram,  in  contrast  to  the  relatively  constant  “background”  (middle 
portion).  The  top  portion  is  due  to  aliasing,  and  will  be  further  discussed, 
below.  Also,  note  the  characteristic  systolic  peak  and  dicrotic  notch.  It  is  im¬ 
portant  to  note  that  the  amplitudes  shown  in  this  figure  are  all  on  logarithm 
scale.  Indeed,  all  of  the  following  analysis  is  performed  on  the  logarithm 
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Figure  2:  The  observed  spectrogram  for  one  patient. 


of  the  observed  spectrogram  amplitude.  The  logarithmic  transformation  is 
motivated  by  the  empirical  result  that  the  logarithm  of  the  amplitudes  is 
bell-shaped. 


3  Methods 


In  this  section  details  of  the  flowchart  (Fig.  1)  are  described. 


3.1  Parameter  Intialization 

The  data  collected  as  a  spectrogrm  are  over  the  domain  of  time  and  doppler 
velocity.  We  develop  a  model  on  the  domain  of  two  independent  variables, 
heart  rate  and  phase  within  a  cardiac  cycle.  The  model  is  posed  this  way 
because  it  is  assumed  that  the  system  of  interest  may  be  approximated  at 
any  instant  in  time  by  a  perfectly  periodic  system.  The  interval  of  idealized 
periodicity  is  called  a  cardiac  cycle.  The  phase,  9,  of  an  event  specifies  the 
time  at  which  that  event  occurs  within  a  cardiac  cycle.  Heart  rate,  u>,  is 
the  number  of  cardiac  cycles  that  would  occur  per  unit  time  in  the  perfectly 
periodic  approximation.  The  time  evolution  of  9  and  cu  are  related  and 
assumed  to  be  stochastic.  The  usual  relation  between  the  two  independent 
variables,  namely  ^  =  u  will  be  posed  as  a  transition  probability  in  a  markov 
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chain  model  of  dynamic  variables  within  the  RBE. 

Each  of  the  components  of  the  proposed  model  introduces  a  number  of 
parameters  (Table  1),  which  are  treated  in  two  different  ways.  The  values  of 
most  of  the  parameters  are  fixed  in  time  for  a  given  patient,  and  are  optimized 
over  the  entire  time  interval  of  the  patient’s  observed  spectrogram.  Two  of 
the  parameters,  however,  are  not  only  dynamic  in  that  they  are  allowed  to 
vary  in  time,  their  possible  values  are  determined  a  priori  over  a  grid.  The 
grid  coordinates  are  the  values  of  the  two  dynamic  parameters  and  are  the 
phase  within  a  cardiac  cycle,  9,  and  heart  rate,  u.  As  a  result,  at  any  time, 
the  aforementioned  pdf  is  a  function  over  these  two  dynamic  parameters. 


3.2  Circulatory  System  Calculation 

A  simple  Windkessel  model,  i.e.  an  electric  circuit  analogy,  is  used  to  produce 
a  time  series  of  total  flow,  or  of  some  other  gross  property  of  flow.  The 
nonlinear  circuit  illustrated  in  Fig.  3.2  exhibits  a  variety  of  realistic  wave 
forms  when  it  is  driven  by  a  square  wave. 

James,  This  is  one  place  where  we  can  use  your  help.  The  bottom  circuit 
diagram  is  what  I  have  been  working  with.  However,  the  results  don’t  look 
too  good.  I  will  continue  working  on  it,  but  because  of  time  constraints 
I  had  to  find  something  that  does  work.  The  top  diagram  seems  to  work 
OK.  So,  the  question  for  you  is  whether  you  can  impose  a  physiological 
interpretation  on  the  top/simple  diagram?  For  the  longer  term,  we  can 
explore  other  diagrams  too,  but  for  now  I’m  stuck  with  just  these  two. 

The  circuit  is  a  rough  analogy  to  a  single-chambered  heart  driving  circu¬ 
lation  through  a  simple  resistive  body.  The  heart  in  this  model  consists  of  a 
lone  ventricle  with  inlet  and  outlet  valves,  and  an  aorta.  The  body  consists 
of  a  single  resistor  representing  all  arteries,  capillaries,  veins,  etc.  The  ac¬ 
tion  of  the  ventricle  is  represented  as  a  voltage  source  with  a  relatively  high 
output  impedance.  The  valves  are  modeled  as  Shockley  diodes  with  finite 
resistance  and  compliance,  representing  the  drag  as  blood  moves  through  the 
valves  and  the  oscillation  of  the  valve  membrane.  The  aorta  is  assumed  to 
be  elastic,  and  the  blood  flowing  into  it  is  assumed  to  have  some  momentum. 
A  simple  resistor  represents  the  flow  through  the  body. 
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The  presence  of  the  diodes  is  essential  in  assuring  realistic  flow  wave  forms. 
Without  the  diodes  in  Fig.  3.2  the  problem  is  linear,  and  so,  can  be  solved  by 
computing  a  transfer  function  mapping  the  driving  square  wave  to  the  cor¬ 
responding  current  through  the  body  resistor,  but  important  features  would 
be  missing  from  the  resulting  wave  form.  With  diodes  the  current  through 
the  resistor  is  dictated  by  a  system  of  first  order  differential  equations,  which 
can  be  solved  numerically. 

The  twelve  parameters  associated  with  the  various  components  in  the  cir¬ 
cuit  analogy  are  tabulated  in  Table  1.  Two  of  these  parameters  -  heart  rate 
(u),  and  phase  in  cardiac  cycle  (0)  -  are  dynamic,  while  the  remaining  ten 
are  fixed. 


parameter 

symbol 

notes  about  prior 

heart  rate 

CO 

gridded 

phase  in  cardiac  cycle 

6 

gridded 

diode  saturation  current  s 

lognorm 

diode  exponent  coeff 

k 

lognorm 

capacitance  of  valve 

cv 

lognorm 

resistance  of  valve 

rv 

lognorm 

inductance  of  aorta 

la 

lognorm 

capacitance  of  aorta 

Ca 

lognorm 

source  impedence 

Th 

lognorm 

source  voltage 

VQ 

lognorm 

square  wave  width 

dwell 

lognorm 

resistance  of  body 

n 

lognorm 

Table  1:  The  dynamic  and  fixed  parameters  introduced  in  the  Circulatory 
System  Calculation  component  of  the  model. 

For  each  gridded  value  of  the  dynamic  parameters  the  circuit  simulation  is 
run  for  a  sufficiently  long  time  so  as  to  remove  transient  behavior.  As  such, 
the  entire  circuit  model  may  be  represented  as  a  gridded  function  f(u,9), 
where  each  gridded  value  corresponds  to  a  maximum  flow  velocity. 

Fig.  3.2  shows  the  time  series  of  maximum  flow  velocity  for  the  patient 
whose  spectroram  is  shown  in  2.  These  time  series  has  some  desirable  fea¬ 
tures,  such  as  the  presence  of  a  dicrotic  notch,  the  realistic  relative  scale 
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Figure  4:  An  example  of  the  time  series  of  maximum  flow  velocity  produced 
by  the  circuit  shown  in  Figure  3.2 


between  the  dicrotic  notch  and  the  systolic  peak,  and  the  qualitative  “expo¬ 
nential”  decay  following  the  peak.  Pierre,  James,  do  you  want  to  say  more 
about  this  figure?? 


3.3  Model  of  Flow  in  the  Middle  Cerebral  Artery:  Flow 
Field  Calculation 

The  aforementioned  windkessel  model  produces  a  time  series  of  total  flow 
through  the  artery.  This  1-D  time  series  is  then  used  as  a  boundary  condition 
which  the  flow  velocity  field  solutions  must  satisfy  in  the  flow  field  model. 

In  the  flow  field  calculation  the  middle  cerebral  artery  is  modeled  as  an 
elastic  tube  ([3]  and  [5]).  The  tube  is  assumed  to  be  straight;  this  assumptions 
is  reasonable  because  although  the  middle  cerebral  artery  is  convoluted,  the 
typical  radius  of  curvature  along  the  artery  is  larger  than  the  radius  of  the 
artery  itself.  The  vessel  wall  is  assumed  to  be  thin,  allowing  one  to  treat  it 
as  a  Hookean  membrane[l]. 

The  mathematical  details  of  this  model  are  similar  to  that  developed  in 
Chapter  17  of  [1].  Specifically,  beginning  with  Equations  17.128,  17.129,  and 
17.131  it  is  straightforward  to  show 
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0  —  Qak2  +  Qbk  +  Qc, 


1  —  v/2  (  —iEh\ 
l 3u>  \cj2paQ/ 


Eh  —i 
2 pa0u>2  oj 


l  —  v/2  iv  —Eh  2  /3p 

/3u>  ao  2paQU! 2  Eh 


Qc  =  (v2 


i'i  &E 
'  Eh 


(1) 

(2) 

(3) 

(4) 


where  ??  make  sure  all  the  symbols  are  defined.??  The  reason  these  questions 
are  derived  here  is  that  equations  17.134  and  17.135  in  [1]  are  erroneous.1 


The  vessel  wall  and  fluid  velocity  solutions  are  required  to  satisfy  the  fol¬ 
lowing  conditions.  The  flow  through  the  artery  is  assumed  to  be  pulsatile, 
and  we  seek  solutions  for  radial  and  axial  velocities  of  the  form  v  (r,  t)  and 
u  (r,  t ),  where  r  is  the  distance  from  the  axis.  As  such,  the  solutions  are  rota- 
tionally  symmetric,  and  there  is  no  angular  velocity.  Furthermore,  the  fluid 
is  assumed  to  be  incompressible  and  satisfying  the  Navier-Stokes  equation. 


Flow  fields  like  this  have  been  studied  by  [3]  and  [5],  and  solutions  are  a 
current  area  of  research,  for  example  in  [6]  and  [15].  Solutions  are  generally 
obtained  through  Fourier  decomposition  of  the  forcing.  In  order  to  arrive  at 
analytic  solutions,  certain  approximations  are  made,  which  in  turn  allow  for 
two  modes  (i.e.  types  of  solutions)  to  exist  for  each  harmonic  [1].  The  most 
general  solution,  a  combination  of  these  two  modes,  therefore  requires  two 
parameters  to  be  specified  uniquely.  The  two  parameters  are  taken  to  be 
the  phase  and  amplitude  at  each  frequency.  The  relative  phase  and  relative 
amplitude  of  the  two  modes  are  assumed  to  be  constant  in  frequency  and 
time.  These  values  are  constrained  to  assure  an  agreement  between  the 
solution  and  the  1-D  time  series  produced  by  the  circulatory  system  model. 
For  a  frequency  zero  (steady  state)  component  it  is  assumed  that  the  axial 
flow  is  a  Hagen-Poiseuille  (parabolic)  type,  and  there  is  no  radial  flow.  2 

1  ??  Do  we  need  to  derive  this?? 

2  extending  the  pulsitile  solution  a  steady  state  solution  would  be  inconsistent  with  the 
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The  basic  physical  equations  defining  the  flow  field  model  [1]  involve  the 
fixed  parameters  in  table  2. 

parameter  symbol  notes  about  prior 


density  of  blood  p 

viscosity  of  blood  p 

density  of  vessel  wall  pw 

equilibrium  radius  of  vessel  a 

thickness  of  vessel  wall  h 

Possion’s  ratio  of  vessel  wall  vw 

Young’s  modulus  of  vessel  wall  E 

relative  amplitude  of  ”  +”  mode  X+ 

relative  amplitude  of  mode  X_ 


? 

? 

? 

? 

? 

? 

? 

the  complex  ratio  between  the  two  modes 
might  be  considered  four  real  numbers 


Table  2:  The  fixed  parameters  introduced  in  the  Flow  Field  Calculation 
component  of  the  model. 


As  seen  in  Fig.  3.3,  the  lowest  Womersley  number  (top  row)  modes  resem¬ 
ble  a  combination  of  Hagen-Poiseuille  flow  and  solid  vibration.  Note  that  as 
the  frequency  increases,  the  shape  of  the  axial  flow  velocity  becomes  increas¬ 
ingly  complex  for  the  “minus”  mode  (top  panel),  and  the  radial  component 
grows  in  the  “plus”  mode  (bottom  panel). 

The  output  of  model  of  blood  flow  in  the  middle  cerebral  artery  is  shown 
in  Fig.  3.3.  The  axes  are  radius  and  time,  with  the  axial  velocity  shown 
parallel  to  the  time  axis.  Typically,  the  radial  velocity  is  on  the  order  of 
fOO  to  200  times  less  than  the  axial  velocity;  here,  they  have  been  scaled 
to  achive  visual  clarity.  It  may  be  tempting  to  ignore  the  lower  amplitude 
radial  component  in  favor  of  the  higher  amplitude  axial  component;  however, 
because  the  observed  spectrograms  span  the  full  domain  of  velocity,  including 
near-zero  velocity,  the  scale  disparity  between  the  two  components  does  not 
justify  discarding  the  lower  amplitude  velocities. 

premise  of  that  flow  model.  In  this  simple  iteration  nothing  is  holding  the  vessel  in  place 
other  than  symmetric  forcing. 
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Figure  5:  An  example  of  the  two  degenerate  modes  of  wave  propogation. 
The  top  (bottom)  panel  corresponds  to  the  “minus”  (’’plus”)  mode,  and  the 
three  rows  in  each  panel  correspond  to  three  values  of  Womersley  number  in 
order  of  increasing  frequency  from  top  to  bottom. 

3.4  Model  of  Transcranial  Doppler  Ultrasound:  Spec¬ 
trogram  Calculation 

The  preceding  steps  in  the  proposed  model  produce  a  2-D  flow  velocity  field. 
TDU,  however,  produces  a  spectrogram,  i.e.,  a  time  series  of  the  histogram  of 
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0.4 


Figure  6:  The  output  of  model  of  blood  flow  in  the  middle  cerebral  artery. 

the  volume  corresponding  to  a  doppler  velocity  v  (where  volume  refers  to  the 
total  number  scatterers  [10]).  In  order  to  compare  this  simulated  blood  flow 
velocity  field  to  the  observed  doppler  ultrasound,  therefore,  a  spectrogram 
must  be  derived  from  the  modeled  velocity  field. 

The  TDU  considered  here  interrogates  only  the  flow  velocity  field  com¬ 
ponent  along  the  transducer  axis.  Though  the  simulated  blood  flow  is  axi- 
symmetric  the  transducer  is  not  assumed  to  be  aligned  with  the  middle  cere¬ 
bral  artery.  Introducing  an  off  bore  angle  $,  the  speed  toward  the  transducer 
v'  is  given  by 

v'  (r,  0,  <b)  =  vaxiai  (r)  cos  <f>  +  vradiai  (r)  sin  <F  cos  0  ,  (5) 

with  0  the  angle  in  cross  section  of  the  artery.  Said  differently,  in  a  cylindri¬ 
cal  coordinate  system  aligned  with  the  artery  r,  0  and  z,  0  is  the  polar  angle 
(around  the  symmetry  axis  of  the  artery),  and  $  is  difference  in  angle  be¬ 
tween  the  symmetry  axis  of  the  ultrasond  foucs  and  the  artery.  This  doppler 
velocity  v'  is  a  function  of  the  corridnate  system  of  artery.  An  ideal  spec¬ 
trogram  will  some  measure  this  v'  but  the  spectrogram  itself  is  a  function  of 
v.  An  ideal  spectrogram  is  a  measure  of  how  much  volume  within  the  focus 
area  is  moving  at  a  velocity  v;  i.e.,  for  what  area  in  the  coordinate  system 
of  the  artery  is  v'  =  v.  Denoting  (r,  0)  with  x ,  and  supressing  <E>,  an  ideal 
spectrogram  s^eai  (v)  can  be  defined  as 

Sideai  (v)  =  lim  —  Area  in  x  such  that  \v'  (x)  —  v\  <  e  (6) 

e— >o  2e 
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Said  differently,  the  ideal  spectrogram  is  the  infinitesimal  area  for  which  the 
doppler  component  of  the  flow  velocity  field  is  equal  to  the  ideal  spectrogram 
velocity.  It  is  relatively  straightforward  to  show  that  for  a  broad  class  of  v', 
the  limit  in  equation  6  can  be  written  as 

s”“ (t,) = L  fwi  <7) 

Where  l  is  any  parameterization  of  contours  for  which  v'  (x)  =  v.  The  re¬ 
lationship  between  an  ideal  spectrogram  and  a  flow  velocity  field  has  been 
studied  by  Evans  [9],  but  not  in  integral/closed  form  as  in  equation  7. 

Additionally,  the  ideal  spectrogram  and  the  observed  spectrogram  are  as¬ 
sumed  to  be  related  by  the  following  statistical  model 

Sdata(v,0(t),uj(t))  =  Sideai(v,6(t),Uj(t))  fv(v)  +  B0,  (8) 

z(v,t)  =  Sdata  (' V,9(t),u(t))  +  N  (v,t )  (9) 


where  A  is  a  zero  mean  gaussian  random  variable,  and  is  a  fixed  parameter 
representing  background  clutter.  The  observed  spectrogram  data  is  z,  and  it 
differs  from  the  rescaled  spectrogram  sdata  by  the  addition  of  noise  N .  The 
variance  of  N  for  a  given  time  and  spectrogram  velocity  is  chosen  a  priori 
during  the  developement  of  this  statistical  model,  the  covariance  of  N  at 
differing  times  and  places  is  assumed  to  be  zero.  An  initial  and  a  refined 


estimate  of  Bo  is  made  for  each  outer  loop  completed  in  Fig.  1. 

parameter 

symbol 

notes  about  prior 

off  bore  angle  of  transducer 

$ 

2d  norm  rayleigh?  30°  thing? 

reflectivity  of  blood 

Bb 

MLE  so  matters  in  tracker  not  in  feedback 

reflectivity  of  tissue 

Bt 

MLE  so  matters  in  tracker  not  in  feedback 

background  clutter 

Bo 

MLE  so  matters  in  tracker  not  in  feedback 

variance  of  noise 

a 

MLE  so  never  matters 

Table  3:  parameters  introduced  in  the 

Spectrogram  Calculation  component 

of  the  model. 
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3.5  Recurssive  Bayesian  Estimation:  Comparison 


The  preceding  analysis  leads  to  an  ideal  spectrogram  depending  on  the  value 
of  the  dynamic  parameters.  Within  a  recursive  Bayesian  framework  the 
observed  spectrogram  z  is  combined  with  the  predicted  spectrogram  Sdata 
to  produce  a  posterior  PDF  over  the  dynamic  parameters  that  evolves  with 
time.  The  location  of  the  maximum  value  of  this  posterior  PDF  is  taken 
as  an  estimate  of  the  values  of  the  dynamic  parameters  as  they  change  over 
time. 


The  dynamic  parameters  and  the  spectrogram  are  discretized  as  follows. 
First,  the  phase  in  the  cardiac  cycle  and  heart  rate  are  both  divided  into  i 
time  steps:  #*,  #j_i, ...,  Oq  and  a;,,  ...,  wo,  respectively.  The  observed  spec¬ 

trogram  is  a  time  series  i),  zC- 1,  •••,  io,  where  z  {j,  i )  is  the  observed  amplitude 
at  the  jth  doppler  velocity  and  ith  time. 

All  fixed  parameters  are  eponymously  fixed  and  a  hypothetical  joint  dis¬ 
tribution  over  all  past  and  present  values  of  the  dynamic  parameters  and 
the  observed  data  is  considered.  The  entire  state  of  the  system  can  then  be 
described  by  a  joint  PDF  of  the  form 

•  (10) 

Given  that  a  spectrogram  is  observed,  the  relevant  quantity  is  the  conditional 
PDF 

P(6i,ei-1,...,90,uji,uJi-1,...,oj0\zi,Zi11,...,z0)  .  (11) 

This  PDF  is  a  function  of  past  and  present  values  of  the  dynamic  parameters, 
but  the  more  useful  quantity  is  the  PDF  of  only  the  latter.  Such  a  PDF  can 
be  obtained  by  marginalizing  over  all  past  values  of  the  dynamic  parameters: 


P  {Oii  A— i)  •••,  At)  (12) 

f  P  ( Oi ,  1,  ...,  Oq,  Ldi,Ldi_  1,  ...,U>o\Zi,  Zi-I,  ...,  Zq)  •  •  -du)^l$) 


which  using  Bayes’  theorem  can  be  rewritten  as 


P  ( Zi , . . .  |  Oi , . . . ,  u>i ,  •  •  • )  P  {Oi , . . . ,  u>i , . . . ) 
P 


dOi-i.-.dOo,  dui-i,  ...du>o 


(14) 


The  denominator  in  Eq.  14  contains  no  parameters,  and  so,  can  be  ignored. 
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In  order  to  estimate  these  PDFs,  statistical  models  are  used.  Specifically, 
the  first  term  in  the  integrand’s  numerator  in  Eq.  14  (often  called  the  likeli¬ 
hood)  is  assumed  to  be 

i 

P  =  JJP(4/|6>j/,uv)  (15) 

i'=  0 

where  each  of  the  PDFs  in  the  product  is  assumed  to  be  multivariate  normal 
with  zero  means  and  identical  covariances: 

P  (zV  1 Ojj ,  (jJj/ )  N  (Zi>  Sdata  (6(/ ,  LUji ) )  ,  (16) 

where  the  components  of  the  vector  $(Mo.  are  given  by  Eq.  8.  The  product  of 
univariate  normal  PDFs  in  Eq.  15  reflects  the  conditional  independence  of 
the  observed  spectrogram  across  time.  The  use  of  if  rather  than  i  indicates 
that  this  relation  holds  for  all  values  of  the  time  index. 

The  second  term  in  the  integrand’s  numerator  in  Eq.  14  (often  called  the 
prior)  involves  only  the  parameters,  and  can  be  written  as  the  product  of 
conditional  PDFs  involving  the  value  of  the  parameters  at  an  earlier  time 
step: 

P  (i Oi ,  ...)  =  P  ( Oi ,  U>i-l)  P  -P  ($0>  ^d-X) 

i 

=  P  (#0>  wo)  JJ  P  ifii'i  kp-l)  (18) 

i'= 1 

The  term  P  (90,u>0)  is  set  to  a  constant. 

Each  of  the  PDFs  in  the  product  Eq.  18  is  assumed  to  be  of  the  form 

P  (0^/ ,  | i ,  LOi' — i)  N  (uty  (Oi'  1  ^if-l))  i  (16) 

where  N  and  V  denote  the  Normal  and  Von  Mises  distributions,  respectively. 
The  argument  of  the  latter  is  a  reflection  of  the  relationship  between  the 
dynamic  variables  lu  ~ 
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Substituting  equations  16  and  19  into  15,  18  and  then  into  14, 


P  (6 UJi  |  Z\—  i ,  . . . ,  Zq)  ^  |  P  (Zi*  |  i  k-V)  1 1  P  ->  |^i'  —  1 5  —  1 )  j  dOi- 

^  Vi^O  i'=l  / 

(20) 

P  ^i)  I  |  |  P  (Zi*  \  @i*  i  k-V)  1 1  P  5  |^i'  —  1 5  ^V  —  l)  j  dOi- 

^  Vi^O  i'=l  / 

(21) 

P  \@ij  ^i)  f  P  {Zi'  —  1  |^i'— 1  >  ^i'— 1 )  P  {@i' ,  (^V  |^i'— 1  ,  ^i'— 1  d9{- 


\i'=l 


(22) 


I  P  {Oh  \@i— 1  5  l)  P  [Zi— 1  |^z— 1  ^  l)  (  I  . .  .ddi—2d^i—2 

Ji- 1  Wz-2 

(23) 


The  final  expression  in  Eq  23,  with  the  nested  integrals  depending  only  on 
the  adjacent  values  i,  is  why  this  method  is  called  recurssive  estimation. 

The  left  panels  in  Fig.  3.5  show  the  probability  of  data,  given  the  dynamic 
parameters  (i.e.,  the  likehood),  and  the  right  panels  display  the  probability  of 
the  dynamic  parameters,  given  data  (i.e.,  the  posterior  PDF).  Note  that  the 
likelihood  does  not  appear  to  have  a  unique  maximum.  In  particular,  there 
is  larger  uncertainty  in  u  than  in  9.  By  contrast,  the  posterior  PDF  for  one 
patient  (top/right  panel)  has  a  unique  maximum,  and  the  uncertainty  is  more 
in  the  9  than  in  cu.  The  posterior  PDF  for  the  other  patient  (bottom/right) 
displays  at  least  two  modes.  Do  we  want  to  explain  this??  The  slanted 
feature  seen  in  the  posterior  PDF  is  characteristic  to  this  method,  because 
it  is  a  consequence  of  the  expression  ^  =  co  discussed  previously. 

In  summary,  the  forward  model  produces  an  ideal  spectrogram  at  gridded 
values  of  9  and  u>.  The  RBE  approach,  then,  allows  one  to  compute  the  con¬ 
ditional  PDF  of  these  parameters,  given  an  observed  spectrogram.  Finally, 
the  dynamic  parameters  are  estimated  to  be  those  values  which  maximize 
that  PDF. 


d9i-iduji-i 
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Figure  7:  The  likelihood  (left)  and  the  posterior  PDF  (right)  over  u,0,  for 
two  patients  (top  and  bottom). 

3.6  Feature  Extraction 

The  dynamic  parameters  which  maximize  the  aforementioned  posterior  PDF 
constitute  a  time  series.  Fig.  3.6  shows  the  time  series  for  two  dynamic 
parameters  oj  and  0  as  estimated  according  to  the  mode  of  the  posterior 
PDF.  As  seen,  the  initial  estimates  change  rapidly;  after  a  few  cardiac  cycles 
the  6  paremeter  settles  on  a  cyclic  pattern,  while  the  u  parameter  converges 
onto  a  fixed  value. 

Armed  with  estimates  of  the  dynamic  parameters,  one  can  generate  a  va¬ 
riety  of  useful  quantities,  such  as  a  spectrogram,  an  envelope,  mean  flow, 
and  even  pulsutility  index.  An  example  of  a  predicted  spectrogram  is  shown 
in  Figure  3.6,  where  it  is  compared  with  the  corresponding  observed  spec¬ 
trogram.  Evidently,  the  fitted  spectrogram  closely  matches  the  observed 
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Figure  8:  The  output  of  the  RBF. 

spectrogram  in  multiple  ways,  For  example,  ??  compare  the  two  figs  -  di¬ 
astolic  notch  -  body  of  the  spectrogram  -  etc.  including  the  aliasing,  etc. 
?? 


Figure  9:  An  example  of  an  observed  (left)  and  a  predicted  (right)  spectro¬ 
gram. 

A  pulsutility  index  can  be  constructed  as  follows:  The  pulsitility  is  defined 
[16]  over  a  cardiac  cycle  as  the  difference  between  the  systolic  and  diastolic 
flow  velocity,  divided  by  the  mean  flow  velocity,  where  systolic  and  dias¬ 
tolic  are  synonomus  with  the  maximum  and  minimum  of  flow  velocity.  The 
c o  which  maximizes  the  posterior  PDF  corresponds  to  a  cardiac  cycle  over 
which  the  pulsitility  may  be  computed,  thus  producing  a  time  series  of  in¬ 
stantaneous  pulsitility.  Similarly  any  qunatity  that  may  be  derived  from  the 
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forward  model  may  be  presented  as  a  time  series  corresponding  to  the  most 
likely  dynamic  parameters. 


4  Results 


This  section  applies  the  above  methodology  to  several  patients,  and  diagnoses 
the  spectrograms  in  terms  of  the  more  fundamental  quantities  such  as  flow 
fields,  etc. 

We  should  agree  on  how  to  interpret  all  of  these  results.  So,  for  now,  I’ll 
say  nothing?? 

I  will  also  attempt  to  produce  some  figures  that  compare  the  predicted  and 
observed  spectrogram  in  a  more  quantitative  fashion.  This  comparison  turns 
out  to  be  more  complicated  than  might  seem?? 

Although  one  aim  of  the  present  work  has  been  to  generate  a  predicted 
spectrogram,  the  flow  field  solution  is  the  more  desirable  (or  fundamental) 
quantity.  However,  there  exists  a  degeneracy  in  the  solutions  that  potentially 
prevents  one  from  obtaining  a  unique,  physical  solution.  The  degeneracy  can 
be  traced  to  an  approximation  made  to  assure  that  a  time  wise  harmonic 
decomposition  of  the  boundary  conditions  and  solutions  is  appropriate  ([1] 
[3],  [5],  [6]).  Each  frequency  component  is  solved  separately,  and  the  final 
solution  is  their  linear  combination  across  frequencies.  In  general,  at  each 
frequency  there  exists  multiple  solutions  (i.e.  degenerate  modes)  of  flow. 
Often  some  of  the  modes  are  dismissed  on  the  grounds  that  they  do  not  meet 
specific  requirements  of  the  problem  at  hand.  Combining  flow  field  solutions 
to  match  a  one  dimensional  time  series  boundary  condition  is  possible  with 
a  linear  combination  of  the  degenerate  modes  for  any  frequency.  Often  in 
the  liturature  an  appeal  to  experiment  is  made  to  resolve  the  true  flow  field 
[17].  In  this  work  (following  [1])  there  are  only  two  modes  present,  and 
the  degeneracy  is  manifested  in  the  relative  phase  and  amplitude  between 
the  two  modes,  at  each  frequency.  An  estimate  of  this  relative  phase  is 
obtained  by  assuming  that  it  is  constant  across  all  freqencies.  As  such,  the 
spectrogram  provides  sufficient  information  to  resolve  the  ambiguity  (because 
it  is  dependent  on  the  entire  velocity  field),  thus  providing  a  unique  flow  field 
solution. 
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Figure  10:  Observed  spectrogram  (left),  flow  field  (middle),  and  predicted 
spectrogram  (right)  for  multiple  patients. 


In  order  to  address  the  sensitivity  of  the  above  results  with  respect  to 
the  choice  of  the  fixed  parameters,  100  different  values  of  all  of  the  fixed 
parameters  were  examined,  and  the  root-mean-squared  (RMS)  error  between 
the  predicted  and  observed  spectrogram  was  computed.  4  shows  three  of  the 
settings  which  produced  relatively  low  values  of  the  RMS  error.  The  columns 
have  been  positioned  in  increasing  order  of  RMS  error  from  left  to  right.  Need 
to  talk  about  these?? 

Overlaid  in  the  top  row  is  the  maximum  value  of  predicted  spectrogram 
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(i.e.,  the  envelope)  in  white.  It  can  be  seen  that  there  are  small  differences  in 
the  three  envelopes.  Indeed,  the  envelope  in  the  middle  panel  appears  to  be 
more  consistent  with  the  underlying  observed  spectrogram.  This  highlights 
the  fact  that  the  producing  spectrograms  which  best  model  the  observed 
spectrogram  may  not  lead  to  better  envelopes.  The  next  section  discusses  an 
alternative  approach  if  the  production  of  an  envelope  (and  not  the  spectro¬ 
gram)  is  the  ultimate  goal. 

As  discussed  previously,  a  compairson  of  the  predicted  and  observed  spec¬ 
trograms  can  aid  in  assessing  the  quality  of  a  flow  field,  even  though  one 
cannot  directly  observe  a  flow  field.  For  example,  as  seen  in  the  right-most 
column  of  Fig.  4,  the  disagreement  between  the  predicted  and  observed 
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Figure  11:  Observed  spectrogram  (top  row),  predicted  spectrogram  (middle 
row),  and  flow  field  (bottom  row)  for  one  patient,  but  with  three  different 
settings  for  the  fixed  parameters. 

spectrograms  implies  that  the  flow  field  (shown  in  bottom  row)  is  unrealis¬ 
tic.  More  ?? 


5  Conclusion  and  Discussion 

A  methodology  has  been  proposed  for  modeling  spectrograms.  The  approach 
involves  a  combination  of  physiologically-based  models  and  statistical  mod¬ 
els.  Each  component  of  the  approach  involves  parameters  which  are  esti¬ 
mated  from  observed  spectrograms,  for  a  given  patient.  With  the  parame¬ 
ters  fitted  to  the  data  from  each  patient,  the  approach  produces  a  predicted 
spectrogram,  which  can  then  be  used  for  generating  summary  measures  such 
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as  envelopes. 

Spectrograms  are  useful  in  assessing  autoregulation,  i.e.,  the  response  of 
smooth  muscle  to  the  C02  concentration  in  blood.  Here,  the  Young’s  mod¬ 
ulus,  Poisson’s  ratio,  and  the  rest  diameter  for  the  middle  cerebral  artery, 
are  parameters  which  model  the  smooth  muscle.  Given  that  the  proposed 
method  estimates  these  parameters,  it  can  therefore  be  considered  a  quan¬ 
tification  of  autoregulation. 

In  the  illustration  here,  the  mode  of  the  PDF  (e.g.,  Figure  ??)  was  used 
to  generate  a  predicted  spectrogram.  At  each  point  in  time,  a  maximum 
flow  velocity  can  be  obtained  from  that  spectrogram,  and  the  resulting  time 
series  is  called  an  envelope.  Given  the  probabilistic  nature  of  the  proposed 
methodology,  one  can  can  additionally  produce  a  PDF  of  the  envelope.  Such 
a  PDF  can  generate  not  only  a  single  time  series,  but  it  can  also  provide  a 
measure  of  uncertainty  of  the  envelope.  This  can  be  done  as  follows:  each 
value  of  u  and  6  corresponds  to  a  maximum  flow  velocity  (produced  by  the 
forward  model).  Given  a  PDF  of  u  and  9,  one  can  then  compute  a  PDF 
for  maximum  flow  velocity.  In  this  way,  attention  is  drwan  away  from  the 
specific  values  of  u>  and  9  which  maximize  their  PDF,  and  therefore,  one  can 
compute  a  PDF  for  any  quantity  produced  by  the  forward  model,  e.g.,  total 
flow,  or  pulsitility  index. 

Clinically,  an  observed  spectrogram  is  not  taken  at  face  value,  and  is  sub¬ 
jected  to  an  internal  “mental  model”  that  exists  only  in  a  clinician’s  mind. 
The  model  developed  here  is  a  quantitification  of  that  mental  model.  As 
such,  an  agreement  between  the  observed  and  the  predicted  spectrogram 
provides  the  clinician  an  endoresment  of  the  model  and  its  parameter  values. 
For  example,  in  cases  of  anticipated  vasospasm  the  mental  model  allows  a 
clinician  to  infer  changes  in  the  radius  of  artery,  max  velocity,  and  total  flow 
volume,  all  from  the  observed  spectrogram.  The  predicted  spectrogram  pro¬ 
duced  by  the  proposed  approach  is  implicitly  based  on  these  parameters  and 
their  values,  and  so  is  more  reliable  than  one  produced  by  a  purely  statistical 
model.  Any  uncertainty  inherent  in  the  PDF  over  the  parameters  is  then 
due  to  model  itself. 

Another  novel  feature  of  the  work  here  is  the  articulation  of  a  closed-form 
relationship  between  a  general  two  or  three  dimensional  flow  field  and  an 
ideal  spectrogram,  namely  Equation  7.  The  utility  of  this  expression  is  that 
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it  allows  one  to  compute  a  spectrogram  from  any  flow  field,  regardless  of 
how  it  is  generated??  Pierre??  For  example,  what??  How  can  we  make  this 
paragraph  longer?? 

Here  the  fixed  parameters  have  not  been  optimized  to  fit  the  data.  Instead 
a  PDF  is  produced  only  for  the  dynamic  parameters,  using  the  measurement 
model  and  recursive  bayesian  filter  approach.  The  outer  loop  showin  in  Fig¬ 
ure  ??  is  a  conceptual  representation  of  the  desire  to  estimate  and  optimize 
the  fixed  parameters.  Iterating  through  the  outerloop  is  a  computationally 
intensive  task.  Specifically,  for  a  single  patient,  one  iteration  requires  about 
1  minute  of  computation  on  a  single-CPU  machine.  Assuming  each  of  the 
parameters  is  sampled  10  times,  a  model  of  the  type  proposed  (i.e.,  with 
about  20  parameters)  will  then  require  1020  minutes,  and  is  therefore  en¬ 
tirely  unfeasible. 

In  this  paper,  for  the  purpose  of  demonstrating  the  approach,  the  outer 
loop  was  not  iterated  a  significant  number  of  times.  In  applying  the  ap¬ 
proach  to  a  single  patient,  it  was  feasible  to  examine  hundreds  of  iterations. 
But  in  applying  it  to  all  of  the  97  patients  in  the  data  set,  only  tens  of 
outer  iterations  were  performed.  A  general  solution  to  this  problem  would 
be  to  optimize  each  component  of  the  model  independently.  For  example,  for 
the  Circulatory  System  Calculation,  the  spectrogram  could  be  preprocessed 
with  and  an  envelope  detection  method  to  produce  a  time  series  of  maximum 
flow  along  with  an  uncertainty  about  that  time  series.  As  this  component 
is  posed  as  a  system  of  first  order  differential  equations  in  terms  fixed  and 
dynamic  paramteters,  it  is  natural  to  employ  a  Kalman  filter.  In  this  al¬ 
ternative  approach  where  the  components  are  independently  optimized,  the 
computational  requirements  would  be  significantly  lower,  i.e.,  of  the  order  of 
10  x  20  minutes  as  compared  to  1020  minutes  noted  above. 

The  modular  nature  of  the  proposed  method  allows  for  customized  predic¬ 
tions.  For  example,  in  a  multi-sensor  situation  where  two  spectrograms  are 
produced  for  the  same  focal  location  but  from  different  directions,  only  the 
“spectrogram  calculation”  component  of  the  model  need  be  changed.  The 
circulatory  system  calculation,  the  flow  field  calculation,  and  the  recursive 
bayesian  filter  comparison  are  unaffected.  The  inputs  to  the  new  components 
would  involve  one  flow  field  from  the  flow  field  calculation  and  two  observed 
spectrograms.  The  output  would  be  the  conditional  probability  of  the  two 
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observed  spectrograms,  given  the  flow  field  and  the  parameters  of  the  new 
component. 

Similarly  the  modular  nature  of  the  model  allows  for  different  flow  velocity 
field  calculations.  For  instance,  closed  form  solutions  to  pulsitile  fluid  flow 
in  an  elastic  tube  date  from  the  50 ’s  in  [5]  and  [3].  But  more  recent  work 
develops  analytic  solutions  for  a  curved  artery  [6].  As  such,  one  can  test 
Myers’  model  by  simply  replacin  the  component  for  flow-field  calculation 
with  Myers’  solutions. 

It  is  worth  mentioning  that  for  the  sake  of  simplicity  it  has  been  assumed 
that  the  focus  volume  of  the  ultrasound  is  a  simple  cross  section  through  the 
artery.  This  assumption  constitutes  an  approximation  to  the  fact  that  the 
ultrasound  focus  volume  has  a  more  complex  spatial  extent  [9] .  One  manner 
in  which  the  assumption  can  be  relaxed  is  by  introducing  a  weighting  function 
in  the  numerator  of  the  integrand  of  Equation  7  to  account  for  the  spatial 
focus  of  the  transducer. 

The  construction  of  a  model  can  also  aid  the  task  of  locating  the  artery 
during  the  initial  set-up  of  the  transducer.  Specifically,  one  can  optimize 
the  expected  distance  to  the  artery  center  in  an  automated  fashion;  and 
this  would  involve  only  the  measurement  model  component.  In  this  case  the 
B  mode  and  M  mode  search  would  be  combined  with  the  larger  inversion 
problem. 

The  model  proposed  here  includes  tissue  motion,  and  therefore,  allows  one 
to  isolate  the  contribution  (to  the  spectrogram)  of  tissue  motion  from  blood 
flow.  A  consequence  of  this  ability  is  that  one  can  ’’explain”  what  is  of¬ 
ten  referred  to  as  aliasing,  namely  features  at  the  top  of  the  spectrogram 
which  appear  to  reflect  gross  features  at  the  bottom  of  the  spectrogram. 
The  motion  of  the  tissue  produces  negative  flow  velocities,  which  are  tradi¬ 
tionally  ’’wrapped  around”  and  displayed  at  the  top  of  the  spectrogram.  In 
the  proposed  model,  these  negative  velocities  can  be  isolated  and  therefore 
filtered-out,  leading  to  physiologically  consistent  spectrograms.  Figure  ?? 
shows  an  illustration  of  this  utility.  The  middle  figure  qualitativiely  resem¬ 
bles  the  spectrogram  (left  figure) ;  and  the  right  figure  shows  the  spectrogram 
with  the  negative  velocities  filtered  out.  It  is  important  to  emphasize  that 
this  filtering  is  physiologically  based. 
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Abstract 

OBJECTIVE: 

To  investigate  the  relation  between  placental  embolisation  and  the  dias¬ 
tolic  notch  in  the  uterine  artery  flow  velocity  waveform  of  pregnant  ewes 
under  general  anaesthesia. 

METHODS: 

Seven  pregnant  ewes  at  a  gestation  16  to  17  weeks  were  anaesthesized 
and  microbeads  of  gelfoam  were  injected  into  the  uterine  artery;  changes 
in  the  uterine  circulation  were  assessed  by  Doppler  velocimetry. 
RESULTS: 

Gelfoam  embolisation  reduced  uterine  blood  flow  in  a  dose-dependent 
manner,  from  a  mean  (95%,  Cl)  of  568  mL/min  (495-641)  to  159  mL/min 
(131-187)  after  the  injection  of  30  mg  of  gelfoam,  and  increased  the  uter¬ 
ine  vascular  resistance  from  135  mmHg  x  min  x  L(-l) (103-167)  to  498 
mmHg  x  min  x  L(-l)  (422-574).  A  diastolic  notch  in  uterine  artery  flow 
velocity  waveform  was  observed  after  20  mg  to  25  mg  of  gelfoam  in  two 
ewes  and  after  injection  of  30  mg  of  gelfoam  in  all  seven  animals.  Injec¬ 
tion  of  30  mg  of  gelfoam  increased  the  pulsatility  index  to  2.4  (1.9-2. 9) 
from  0.6  (0.5-0. 7).  The  mean  uterine  vascular  resistance  at  the  time  of  the 
appearance  of  a  diastolic  notch  was  414  mmHg  x  min  x  L(-l) (377-451). 
CONCLUSION: 

These  findings  suggest  that  an  elevated  pulsatility  index  and  the  pres- 
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